Histopathology Research Template

serdarbalci.com

---

https://doi.org/10.5281/zenodo.3635430

https://osf.io/3tjfk/

Histopathology Research Template 🔬

---

1 Introduction

• State the marker of interest, the study objectives, and hypotheses (Knijn, Simmer, and Nagtegaal 2015).¹

2 Materials & Methods

Describe Materials and Methods as highlighted in (Knijn, Simmer, and Nagtegaal 2015).²

• Describe patient characteristics, and inclusion and exclusion criteria

• Describe treatment details

• Describe the type of material used

• Specify how expression of the biomarker was assessed

• Describe the number of independent (blinded) scorers and how they scored

• State the method of case selection, study design, origin of the cases, and time frame

• Describe the end of the follow-up period and median follow-up time

• Define all clinical endpoints examined

• Specify all applied statistical methods

• Describe how interactions with other clinical/pathological factors were analyzed

---

2.1 Header Codes

Codes for general settings.³

Setup global chunk settings⁴

knitr::opts_chunk$set(
    eval = TRUE,
    echo = TRUE,
    fig.path = here::here("figs/"),
    message = FALSE,
    warning = FALSE,
    error = FALSE,
    cache = TRUE,
    comment = NA,
    tidy = TRUE,
    fig.width = 6,
    fig.height = 4
)

Load Library

see R/loadLibrary.R for the libraries loaded.

source(file = here::here("R", "loadLibrary.R"))

---

2.2 Generate Fake Data

Codes for generating fake data.⁵

Generate Fake Data

This code generates a fake histopathological data. Some sources for fake data generation here⁶ , here⁷ , here⁸ , here⁹ , here¹⁰ , here¹¹ , here¹² , here¹³ , and here¹⁴ .

Use this code to generate fake clinicopathologic data

source(file = here::here("R", "gc_fake_data.R"))

wakefield::table_heat(x = fakedata, palette = "Set1", flip = TRUE, print = TRUE)

---

2.3 Import Data

Codes for importing data.¹⁵

Read the data

library(readxl)
mydata <- readxl::read_excel(here::here("data", "mydata.xlsx"))
# View(mydata) # Use to view data after importing

Add code for import multiple data purrr reduce

---

2.4 Study Population

2.4.1 Report General Features

Codes for reporting general features.¹⁶

Dataframe Report

# Dataframe report
mydata %>% dplyr::select(-contains("Date")) %>% report::report(.)

The data contains 250 observations of the following variables:
  - ID: 250 entries: 001, n = 1; 002, n = 1; 003, n = 1 and 247 others (0 missing)
  - Name: 249 entries: Aaleyah, n = 1; Abrea, n = 1; Afonso, n = 1 and 246 others (1 missing)
  - Sex: 2 entries: Male, n = 127; Female, n = 122 (1 missing)
  - Age: Mean = 48.59, SD = 14.12, Median = , MAD = 17.79, range: [25, 73], Skewness = 0.07, Kurtosis = -1.17, 1 missing
  - Race: 7 entries: White, n = 151; Hispanic, n = 49; Black, n = 30 and 4 others (1 missing)
  - PreinvasiveComponent: 2 entries: Absent, n = 205; Present, n = 44 (1 missing)
  - LVI: 2 entries: Absent, n = 148; Present, n = 102 (0 missing)
  - PNI: 2 entries: Absent, n = 181; Present, n = 68 (1 missing)
  - Death: 2 levels: FALSE (n = 80, 32.00%); TRUE (n = 169, 67.60%) and missing (n = 1, 0.40%)
  - Group: 2 entries: Treatment, n = 130; Control, n = 119 (1 missing)
  - Grade: 3 entries: 3, n = 97; 2, n = 84; 1, n = 68 (1 missing)
  - TStage: 4 entries: 4, n = 113; 3, n = 65; 2, n = 49 and 1 other (0 missing)
  - AntiX_intensity: Mean = 2.39, SD = 0.66, Median = , MAD = 1.48, range: [1, 3], Skewness = -0.62, Kurtosis = -0.66, 1 missing
  - AntiY_intensity: Mean = 1.91, SD = 0.75, Median = , MAD = 1.48, range: [1, 3], Skewness = 0.15, Kurtosis = -1.20, 1 missing
  - LymphNodeMetastasis: 2 entries: Absent, n = 147; Present, n = 102 (1 missing)
  - Valid: 2 levels: FALSE (n = 134, 53.60%); TRUE (n = 115, 46.00%) and missing (n = 1, 0.40%)
  - Smoker: 2 levels: FALSE (n = 132, 52.80%); TRUE (n = 117, 46.80%) and missing (n = 1, 0.40%)
  - Grade_Level: 3 entries: high, n = 91; moderate, n = 83; low, n = 75 (1 missing)
  - DeathTime: 2 entries: Within1Year, n = 149; MoreThan1Year, n = 101 (0 missing)

mydata %>% explore::describe_tbl()

250 observations with 21 variables
17 variables containing missings (NA)
0 variables with no variance

---

2.5 Ethics and IRB

2.5.1 Always Respect Patient Privacy

Always Respect Patient Privacy
- Health Information Privacy¹⁷
- Kişisel Verilerin Korunması¹⁸

---

2.6 Define Variable Types

Codes for defining variable types.¹⁹

print column names as vector

dput(names(mydata))

c("ID", "Name", "Sex", "Age", "Race", "PreinvasiveComponent", 
"LVI", "PNI", "LastFollowUpDate", "Death", "Group", "Grade", 
"TStage", "AntiX_intensity", "AntiY_intensity", "LymphNodeMetastasis", 
"Valid", "Smoker", "Grade_Level", "SurgeryDate", "DeathTime")

2.6.1 Find Key Columns

2.6.1.1 Find ID and key columns to exclude from analysis

vctrs::vec_assert()

dplyr::all_equal()

arsenal::compare()

visdat::vis_compare()

See the code as function in R/find_key.R.

keycolumns <- mydata %>% sapply(., FUN = dataMaid::isKey) %>% tibble::as_tibble() %>% 
    dplyr::select(which(.[1, ] == TRUE)) %>% names()
keycolumns

[1] "ID"   "Name"

2.6.2 Variable Types

Get variable types

mydata %>% dplyr::select(-keycolumns) %>% inspectdf::inspect_types()

# A tibble: 4 x 4
  type             cnt  pcnt col_name  
  <chr>          <int> <dbl> <list>    
1 character         11  57.9 <chr [11]>
2 logical            3  15.8 <chr [3]> 
3 numeric            3  15.8 <chr [3]> 
4 POSIXct POSIXt     2  10.5 <chr [2]> 

mydata %>% dplyr::select(-keycolumns, -contains("Date")) %>% describer::describe() %>% 
    knitr::kable(format = "markdown")

.column_name	.column_class	.column_type	.count_elements	.mean_value	.sd_value	.q0_value	.q25_value	.q50_value	.q75_value	.q100_value
Sex	character	character	250	NA	NA	Female	NA	NA	NA	Male
Age	numeric	double	250	48.586345	14.1156504	25	36	48	61	73
Race	character	character	250	NA	NA	Asian	NA	NA	NA	White
PreinvasiveComponent	character	character	250	NA	NA	Absent	NA	NA	NA	Present
LVI	character	character	250	NA	NA	Absent	NA	NA	NA	Present
PNI	character	character	250	NA	NA	Absent	NA	NA	NA	Present
Death	logical	logical	250	NA	NA	FALSE	NA	NA	NA	TRUE
Group	character	character	250	NA	NA	Control	NA	NA	NA	Treatment
Grade	character	character	250	NA	NA	1	NA	NA	NA	3
TStage	character	character	250	NA	NA	1	NA	NA	NA	4
AntiX_intensity	numeric	double	250	2.385542	0.6629235	1	2	2	3	3
AntiY_intensity	numeric	double	250	1.907631	0.7483194	1	1	2	2	3
LymphNodeMetastasis	character	character	250	NA	NA	Absent	NA	NA	NA	Present
Valid	logical	logical	250	NA	NA	FALSE	NA	NA	NA	TRUE
Smoker	logical	logical	250	NA	NA	FALSE	NA	NA	NA	TRUE
Grade_Level	character	character	250	NA	NA	high	NA	NA	NA	moderate
DeathTime	character	character	250	NA	NA	MoreThan1Year	NA	NA	NA	Within1Year

Plot variable types

mydata %>% dplyr::select(-keycolumns) %>% inspectdf::inspect_types() %>% inspectdf::show_plot()

# https://github.com/ropensci/visdat
# http://visdat.njtierney.com/articles/using_visdat.html
# https://cran.r-project.org/web/packages/visdat/index.html
# http://visdat.njtierney.com/

# visdat::vis_guess(mydata)

visdat::vis_dat(mydata)

mydata %>% explore::explore_tbl()

2.6.3 Define Variable Types

2.6.3.1 Find character variables

characterVariables <- mydata %>% dplyr::select(-keycolumns) %>% inspectdf::inspect_types() %>% 
    dplyr::filter(type == "character") %>% dplyr::select(col_name) %>% dplyr::pull() %>% 
    unlist()

characterVariables

 [1] "Sex"                  "Race"                 "PreinvasiveComponent"
 [4] "LVI"                  "PNI"                  "Group"               
 [7] "Grade"                "TStage"               "LymphNodeMetastasis" 
[10] "Grade_Level"          "DeathTime"           

2.6.3.2 Find categorical variables

categoricalVariables <- mydata %>% dplyr::select(-keycolumns, -contains("Date")) %>% 
    describer::describe() %>% janitor::clean_names() %>% dplyr::filter(column_type == 
    "factor") %>% dplyr::select(column_name) %>% dplyr::pull()

categoricalVariables

character(0)

2.6.3.3 Find continious variables

continiousVariables <- mydata %>% dplyr::select(-keycolumns, -contains("Date")) %>% 
    describer::describe() %>% janitor::clean_names() %>% dplyr::filter(column_type == 
    "numeric" | column_type == "double") %>% dplyr::select(column_name) %>% dplyr::pull()

continiousVariables

[1] "Age"             "AntiX_intensity" "AntiY_intensity"

2.6.3.4 Find numeric variables

numericVariables <- mydata %>% dplyr::select(-keycolumns) %>% inspectdf::inspect_types() %>% 
    dplyr::filter(type == "numeric") %>% dplyr::select(col_name) %>% dplyr::pull() %>% 
    unlist()

numericVariables

[1] "Age"             "AntiX_intensity" "AntiY_intensity"

2.6.3.5 Find integer variables

integerVariables <- mydata %>% dplyr::select(-keycolumns) %>% inspectdf::inspect_types() %>% 
    dplyr::filter(type == "integer") %>% dplyr::select(col_name) %>% dplyr::pull() %>% 
    unlist()

integerVariables

NULL

2.6.3.6 Find list variables

listVariables <- mydata %>% dplyr::select(-keycolumns) %>% inspectdf::inspect_types() %>% 
    dplyr::filter(type == "list") %>% dplyr::select(col_name) %>% dplyr::pull() %>% 
    unlist()
listVariables

NULL

2.6.3.7 Find date variables

is_date <- function(x) inherits(x, c("POSIXct", "POSIXt"))

dateVariables <- names(which(sapply(mydata, FUN = is_date) == TRUE))
dateVariables

[1] "LastFollowUpDate" "SurgeryDate"     

---

2.7 Overview the Data

Codes for overviewing the data.²⁰

2.7.1 View Data

View(mydata)

reactable::reactable(data = mydata, sortable = TRUE, resizable = TRUE, filterable = TRUE, 
    searchable = TRUE, pagination = TRUE, paginationType = "numbers", showPageSizeOptions = TRUE, 
    highlight = TRUE, striped = TRUE, outlined = TRUE, compact = TRUE, wrap = FALSE, 
    showSortIcon = TRUE, showSortable = TRUE)

2.7.2 Overview / Exploratory Data Analysis (EDA)

Summary of Data via summarytools 📦

summarytools::view(summarytools::dfSummary(mydata %>% dplyr::select(-keycolumns)))

if (!dir.exists(here::here("out"))) {
    dir.create(here::here("out"))
}

summarytools::view(x = summarytools::dfSummary(mydata %>% dplyr::select(-keycolumns)), 
    file = here::here("out", "mydata_summary.html"))

Summary via dataMaid 📦

if (!dir.exists(here::here("out"))) {
    dir.create(here::here("out"))
}

dataMaid::makeDataReport(data = mydata, file = here::here("out", "dataMaid_mydata.Rmd"), 
    replace = TRUE, openResult = FALSE, render = FALSE, quiet = TRUE)

Summary via explore 📦

if (!dir.exists(here::here("out"))) {
    dir.create(here::here("out"))
}

mydata %>% dplyr::select(-dateVariables) %>% explore::report(output_file = "mydata_report.html", 
    output_dir = here::here("out"))

Glimpse of Data

dplyr::glimpse(mydata %>% dplyr::select(-keycolumns, -dateVariables))

Observations: 250
Variables: 17
$ Sex                  <chr> "Female", "Female", "Female", "Female", "Male", …
$ Age                  <dbl> 54, 39, 32, 58, 38, 50, 36, 70, 43, 45, 65, 27, …
$ Race                 <chr> "White", "White", "White", "Hispanic", "White", …
$ PreinvasiveComponent <chr> "Present", "Absent", "Absent", "Present", "Prese…
$ LVI                  <chr> "Present", "Absent", "Present", "Absent", "Absen…
$ PNI                  <chr> "Absent", "Absent", "Absent", "Present", "Absent…
$ Death                <lgl> TRUE, TRUE, FALSE, TRUE, TRUE, FALSE, FALSE, TRU…
$ Group                <chr> "Treatment", "Treatment", "Control", "Treatment"…
$ Grade                <chr> "1", "3", "3", "3", "2", "1", "2", "2", "1", "3"…
$ TStage               <chr> "4", "4", "2", "1", "1", "4", "4", "4", "4", "1"…
$ AntiX_intensity      <dbl> 2, 3, 3, 3, 2, 2, 1, 3, 2, 2, 3, 2, 3, 1, 3, 3, …
$ AntiY_intensity      <dbl> 1, 1, 1, 1, 3, 3, 2, 3, 2, 2, 2, 1, 1, 2, 2, 2, …
$ LymphNodeMetastasis  <chr> "Present", "Absent", "Absent", "Absent", "Presen…
$ Valid                <lgl> FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, F…
$ Smoker               <lgl> FALSE, FALSE, TRUE, TRUE, FALSE, TRUE, FALSE, TR…
$ Grade_Level          <chr> "low", "low", "moderate", "high", "high", "moder…
$ DeathTime            <chr> "Within1Year", "Within1Year", "Within1Year", "Wi…

mydata %>% explore::describe()

# A tibble: 21 x 8
   variable             type     na na_pct unique   min  mean   max
   <chr>                <chr> <int>  <dbl>  <int> <dbl> <dbl> <dbl>
 1 ID                   chr       0    0      250    NA NA       NA
 2 Name                 chr       1    0.4    250    NA NA       NA
 3 Sex                  chr       1    0.4      3    NA NA       NA
 4 Age                  dbl       1    0.4     50    25 48.6     73
 5 Race                 chr       1    0.4      8    NA NA       NA
 6 PreinvasiveComponent chr       1    0.4      3    NA NA       NA
 7 LVI                  chr       0    0        2    NA NA       NA
 8 PNI                  chr       1    0.4      3    NA NA       NA
 9 LastFollowUpDate     dat       1    0.4     13    NA NA       NA
10 Death                lgl       1    0.4      3     0  0.68     1
# … with 11 more rows

Explore

explore::explore(mydata)

2.7.3 Control Data

Control Data if matching expectations

visdat::vis_expect(data = mydata, expectation = ~.x == -1, show_perc = TRUE)

visdat::vis_expect(mydata, ~.x >= 25)

See missing values

visdat::vis_miss(airquality, cluster = TRUE)

visdat::vis_miss(airquality, sort_miss = TRUE)

xray::anomalies(mydata)

$variables
               Variable   q qNA  pNA qZero pZero qBlank pBlank qInf pInf
1                 Valid 250   1 0.4%   134 53.6%      0      -    0    -
2                Smoker 250   1 0.4%   132 52.8%      0      -    0    -
3                 Death 250   1 0.4%    80   32%      0      -    0    -
4                   Sex 250   1 0.4%     0     -      0      -    0    -
5  PreinvasiveComponent 250   1 0.4%     0     -      0      -    0    -
6                   PNI 250   1 0.4%     0     -      0      -    0    -
7                 Group 250   1 0.4%     0     -      0      -    0    -
8   LymphNodeMetastasis 250   1 0.4%     0     -      0      -    0    -
9                 Grade 250   1 0.4%     0     -      0      -    0    -
10      AntiX_intensity 250   1 0.4%     0     -      0      -    0    -
11      AntiY_intensity 250   1 0.4%     0     -      0      -    0    -
12          Grade_Level 250   1 0.4%     0     -      0      -    0    -
13                 Race 250   1 0.4%     0     -      0      -    0    -
14     LastFollowUpDate 250   1 0.4%     0     -      0      -    0    -
15                  Age 250   1 0.4%     0     -      0      -    0    -
16          SurgeryDate 250   1 0.4%     0     -      0      -    0    -
17                 Name 250   1 0.4%     0     -      0      -    0    -
18                  LVI 250   0    -     0     -      0      -    0    -
19            DeathTime 250   0    -     0     -      0      -    0    -
20               TStage 250   0    -     0     -      0      -    0    -
21                   ID 250   0    -     0     -      0      -    0    -
   qDistinct      type anomalous_percent
1          3   Logical               54%
2          3   Logical             53.2%
3          3   Logical             32.4%
4          3 Character              0.4%
5          3 Character              0.4%
6          3 Character              0.4%
7          3 Character              0.4%
8          3 Character              0.4%
9          4 Character              0.4%
10         4   Numeric              0.4%
11         4   Numeric              0.4%
12         4 Character              0.4%
13         8 Character              0.4%
14        13 Timestamp              0.4%
15        50   Numeric              0.4%
16       229 Timestamp              0.4%
17       250 Character              0.4%
18         2 Character                 -
19         2 Character                 -
20         4 Character                 -
21       250 Character                 -

$problem_variables
 [1] Variable          q                 qNA               pNA              
 [5] qZero             pZero             qBlank            pBlank           
 [9] qInf              pInf              qDistinct         type             
[13] anomalous_percent problems         
<0 rows> (or 0-length row.names)

xray::distributions(mydata)

================================================================================

[1] "Ignoring variable LastFollowUpDate: Unsupported type for visualization."
[1] "Ignoring variable SurgeryDate: Unsupported type for visualization."

         Variable p_1 p_10 p_25 p_50 p_75 p_90 p_99
1 AntiX_intensity   1  1.8    2    2    3    3    3
2 AntiY_intensity   1    1    1    2    2    3    3
3             Age  25 29.8   36   48   61 68.2   73

2.7.4 Explore Data

Summary of Data via DataExplorer 📦

DataExplorer::plot_str(mydata)

DataExplorer::plot_str(mydata, type = "r")

DataExplorer::introduce(mydata)

# A tibble: 1 x 9
   rows columns discrete_columns continuous_colu… all_missing_col…
  <int>   <int>            <int>            <int>            <int>
1   250      21               18                3                0
# … with 4 more variables: total_missing_values <int>, complete_rows <int>,
#   total_observations <int>, memory_usage <dbl>

DataExplorer::plot_intro(mydata)

DataExplorer::plot_missing(mydata)

Drop columns

mydata2 <- DataExplorer::drop_columns(mydata, "TStage")

DataExplorer::plot_bar(mydata)

DataExplorer::plot_bar(mydata, with = "Death")

DataExplorer::plot_histogram(mydata)

---

---

3 Statistical Analysis

Learn these tests as highlighted in (Schmidt et al. 2017).²¹

---

4 Results

Write results as described in (Knijn, Simmer, and Nagtegaal 2015)²²

• Describe the number of patients included in the analysis and reason for dropout

• Report patient/disease characteristics (including the biomarker of interest) with the number of missing values

• Describe the interaction of the biomarker of interest with established prognostic variables

• Include at least 90 % of initial cases included in univariate and multivariate analyses

• Report the estimated effect (relative risk/odds ratio, confidence interval, and p value) in univariate analysis

• Report the estimated effect (hazard rate/odds ratio, confidence interval, and p value) in multivariate analysis

• Report the estimated effects (hazard ratio/odds ratio, confidence interval, and p value) of other prognostic factors included in multivariate analysis

---

4.1 Descriptive Statistics

Codes for Descriptive Statistics.²³

4.1.1 Table One

Report Data properties via report 📦

mydata %>% dplyr::select(-dplyr::contains("Date")) %>% report::report()

The data contains 250 observations of the following variables:
  - ID: 250 entries: 001, n = 1; 002, n = 1; 003, n = 1 and 247 others (0 missing)
  - Name: 249 entries: Aaleyah, n = 1; Abrea, n = 1; Afonso, n = 1 and 246 others (1 missing)
  - Sex: 2 entries: Male, n = 127; Female, n = 122 (1 missing)
  - Age: Mean = 48.59, SD = 14.12, Median = , MAD = 17.79, range: [25, 73], Skewness = 0.07, Kurtosis = -1.17, 1 missing
  - Race: 7 entries: White, n = 151; Hispanic, n = 49; Black, n = 30 and 4 others (1 missing)
  - PreinvasiveComponent: 2 entries: Absent, n = 205; Present, n = 44 (1 missing)
  - LVI: 2 entries: Absent, n = 148; Present, n = 102 (0 missing)
  - PNI: 2 entries: Absent, n = 181; Present, n = 68 (1 missing)
  - Death: 2 levels: FALSE (n = 80, 32.00%); TRUE (n = 169, 67.60%) and missing (n = 1, 0.40%)
  - Group: 2 entries: Treatment, n = 130; Control, n = 119 (1 missing)
  - Grade: 3 entries: 3, n = 97; 2, n = 84; 1, n = 68 (1 missing)
  - TStage: 4 entries: 4, n = 113; 3, n = 65; 2, n = 49 and 1 other (0 missing)
  - AntiX_intensity: Mean = 2.39, SD = 0.66, Median = , MAD = 1.48, range: [1, 3], Skewness = -0.62, Kurtosis = -0.66, 1 missing
  - AntiY_intensity: Mean = 1.91, SD = 0.75, Median = , MAD = 1.48, range: [1, 3], Skewness = 0.15, Kurtosis = -1.20, 1 missing
  - LymphNodeMetastasis: 2 entries: Absent, n = 147; Present, n = 102 (1 missing)
  - Valid: 2 levels: FALSE (n = 134, 53.60%); TRUE (n = 115, 46.00%) and missing (n = 1, 0.40%)
  - Smoker: 2 levels: FALSE (n = 132, 52.80%); TRUE (n = 117, 46.80%) and missing (n = 1, 0.40%)
  - Grade_Level: 3 entries: high, n = 91; moderate, n = 83; low, n = 75 (1 missing)
  - DeathTime: 2 entries: Within1Year, n = 149; MoreThan1Year, n = 101 (0 missing)

Table 1 via arsenal 📦

# cat(names(mydata), sep = " + \n")
library(arsenal)
tab1 <- arsenal::tableby(
  ~ Sex +
    Age +
    Race +
    PreinvasiveComponent +
    LVI +
    PNI +
    Death +
    Group +
    Grade +
    TStage +
    # `Anti-X-intensity` +
    # `Anti-Y-intensity` +
    LymphNodeMetastasis +
    Valid +
    Smoker +
    Grade_Level
  ,
  data = mydata 
)
summary(tab1)

	Overall (N=250)
Sex
N-Miss	1
Female	122 (49.0%)
Male	127 (51.0%)
Age
N-Miss	1
Mean (SD)	48.586 (14.116)
Range	25.000 - 73.000
Race
N-Miss	1
Asian	11 (4.4%)
Bi-Racial	5 (2.0%)
Black	30 (12.0%)
Hispanic	49 (19.7%)
Native	1 (0.4%)
Other	2 (0.8%)
White	151 (60.6%)
PreinvasiveComponent
N-Miss	1
Absent	205 (82.3%)
Present	44 (17.7%)
LVI
Absent	148 (59.2%)
Present	102 (40.8%)
PNI
N-Miss	1
Absent	181 (72.7%)
Present	68 (27.3%)
Death
N-Miss	1
FALSE	80 (32.1%)
TRUE	169 (67.9%)
Group
N-Miss	1
Control	119 (47.8%)
Treatment	130 (52.2%)
Grade
N-Miss	1
1	68 (27.3%)
2	84 (33.7%)
3	97 (39.0%)
TStage
1	23 (9.2%)
2	49 (19.6%)
3	65 (26.0%)
4	113 (45.2%)
LymphNodeMetastasis
N-Miss	1
Absent	147 (59.0%)
Present	102 (41.0%)
Valid
N-Miss	1
FALSE	134 (53.8%)
TRUE	115 (46.2%)
Smoker
N-Miss	1
FALSE	132 (53.0%)
TRUE	117 (47.0%)
Grade_Level
N-Miss	1
high	91 (36.5%)
low	75 (30.1%)
moderate	83 (33.3%)

Table 1 via tableone 📦

library(tableone)
mydata %>% dplyr::select(-keycolumns, -dateVariables) %>% tableone::CreateTableOne(data = .)

                                    
                                     Overall      
  n                                    250        
  Sex = Male (%)                       127 (51.0) 
  Age (mean (SD))                    48.59 (14.12)
  Race (%)                                        
     Asian                              11 ( 4.4) 
     Bi-Racial                           5 ( 2.0) 
     Black                              30 (12.0) 
     Hispanic                           49 (19.7) 
     Native                              1 ( 0.4) 
     Other                               2 ( 0.8) 
     White                             151 (60.6) 
  PreinvasiveComponent = Present (%)    44 (17.7) 
  LVI = Present (%)                    102 (40.8) 
  PNI = Present (%)                     68 (27.3) 
  Death = TRUE (%)                     169 (67.9) 
  Group = Treatment (%)                130 (52.2) 
  Grade (%)                                       
     1                                  68 (27.3) 
     2                                  84 (33.7) 
     3                                  97 (39.0) 
  TStage (%)                                      
     1                                  23 ( 9.2) 
     2                                  49 (19.6) 
     3                                  65 (26.0) 
     4                                 113 (45.2) 
  AntiX_intensity (mean (SD))         2.39 (0.66) 
  AntiY_intensity (mean (SD))         1.91 (0.75) 
  LymphNodeMetastasis = Present (%)    102 (41.0) 
  Valid = TRUE (%)                     115 (46.2) 
  Smoker = TRUE (%)                    117 (47.0) 
  Grade_Level (%)                                 
     high                               91 (36.5) 
     low                                75 (30.1) 
     moderate                           83 (33.3) 
  DeathTime = Within1Year (%)          149 (59.6) 

Descriptive Statistics of Continuous Variables

mydata %>% dplyr::select(continiousVariables, numericVariables, integerVariables) %>% 
    summarytools::descr(., style = "rmarkdown")

print(summarytools::descr(mydata), method = "render", table.classes = "st-small")

mydata %>% summarytools::descr(., stats = "common", transpose = TRUE, headings = FALSE)

mydata %>% summarytools::descr(stats = "common") %>% summarytools::tb()

mydata$Sex %>% summarytools::freq(cumul = FALSE, report.nas = FALSE) %>% summarytools::tb()

mydata %>% explore::describe() %>% dplyr::filter(unique < 5)

# A tibble: 15 x 8
   variable             type     na na_pct unique   min  mean   max
   <chr>                <chr> <int>  <dbl>  <int> <dbl> <dbl> <dbl>
 1 Sex                  chr       1    0.4      3    NA NA       NA
 2 PreinvasiveComponent chr       1    0.4      3    NA NA       NA
 3 LVI                  chr       0    0        2    NA NA       NA
 4 PNI                  chr       1    0.4      3    NA NA       NA
 5 Death                lgl       1    0.4      3     0  0.68     1
 6 Group                chr       1    0.4      3    NA NA       NA
 7 Grade                chr       1    0.4      4    NA NA       NA
 8 TStage               chr       0    0        4    NA NA       NA
 9 AntiX_intensity      dbl       1    0.4      4     1  2.39     3
10 AntiY_intensity      dbl       1    0.4      4     1  1.91     3
11 LymphNodeMetastasis  chr       1    0.4      3    NA NA       NA
12 Valid                lgl       1    0.4      3     0  0.46     1
13 Smoker               lgl       1    0.4      3     0  0.47     1
14 Grade_Level          chr       1    0.4      4    NA NA       NA
15 DeathTime            chr       0    0        2    NA NA       NA

mydata %>% explore::describe() %>% dplyr::filter(na > 0)

# A tibble: 17 x 8
   variable             type     na na_pct unique   min  mean   max
   <chr>                <chr> <int>  <dbl>  <int> <dbl> <dbl> <dbl>
 1 Name                 chr       1    0.4    250    NA NA       NA
 2 Sex                  chr       1    0.4      3    NA NA       NA
 3 Age                  dbl       1    0.4     50    25 48.6     73
 4 Race                 chr       1    0.4      8    NA NA       NA
 5 PreinvasiveComponent chr       1    0.4      3    NA NA       NA
 6 PNI                  chr       1    0.4      3    NA NA       NA
 7 LastFollowUpDate     dat       1    0.4     13    NA NA       NA
 8 Death                lgl       1    0.4      3     0  0.68     1
 9 Group                chr       1    0.4      3    NA NA       NA
10 Grade                chr       1    0.4      4    NA NA       NA
11 AntiX_intensity      dbl       1    0.4      4     1  2.39     3
12 AntiY_intensity      dbl       1    0.4      4     1  1.91     3
13 LymphNodeMetastasis  chr       1    0.4      3    NA NA       NA
14 Valid                lgl       1    0.4      3     0  0.46     1
15 Smoker               lgl       1    0.4      3     0  0.47     1
16 Grade_Level          chr       1    0.4      4    NA NA       NA
17 SurgeryDate          dat       1    0.4    229    NA NA       NA

mydata %>% explore::describe()

# A tibble: 21 x 8
   variable             type     na na_pct unique   min  mean   max
   <chr>                <chr> <int>  <dbl>  <int> <dbl> <dbl> <dbl>
 1 ID                   chr       0    0      250    NA NA       NA
 2 Name                 chr       1    0.4    250    NA NA       NA
 3 Sex                  chr       1    0.4      3    NA NA       NA
 4 Age                  dbl       1    0.4     50    25 48.6     73
 5 Race                 chr       1    0.4      8    NA NA       NA
 6 PreinvasiveComponent chr       1    0.4      3    NA NA       NA
 7 LVI                  chr       0    0        2    NA NA       NA
 8 PNI                  chr       1    0.4      3    NA NA       NA
 9 LastFollowUpDate     dat       1    0.4     13    NA NA       NA
10 Death                lgl       1    0.4      3     0  0.68     1
# … with 11 more rows

4.1.2 Categorical Variables

Use R/gc_desc_cat.R to generate gc_desc_cat.Rmd containing descriptive statistics for categorical variables

source(here::here("R", "gc_desc_cat.R"))

4.1.2.1 Descriptive Statistics Sex

mydata %>% janitor::tabyl(Sex) %>% janitor::adorn_pct_formatting(rounding = "half up", 
    digits = 1) %>% knitr::kable()

Sex	n	percent	valid_percent
Female	122	48.8%	49.0%
Male	127	50.8%	51.0%
NA	1	0.4%	-

4.1.2.2 Descriptive Statistics Race

mydata %>% janitor::tabyl(Race) %>% janitor::adorn_pct_formatting(rounding = "half up", 
    digits = 1) %>% knitr::kable()

Race	n	percent	valid_percent
Asian	11	4.4%	4.4%
Bi-Racial	5	2.0%	2.0%
Black	30	12.0%	12.0%
Hispanic	49	19.6%	19.7%
Native	1	0.4%	0.4%
Other	2	0.8%	0.8%
White	151	60.4%	60.6%
NA	1	0.4%	-

4.1.2.3 Descriptive Statistics PreinvasiveComponent

mydata %>% janitor::tabyl(PreinvasiveComponent) %>% janitor::adorn_pct_formatting(rounding = "half up", 
    digits = 1) %>% knitr::kable()

PreinvasiveComponent	n	percent	valid_percent
Absent	205	82.0%	82.3%
Present	44	17.6%	17.7%
NA	1	0.4%	-

4.1.2.4 Descriptive Statistics LVI

mydata %>% janitor::tabyl(LVI) %>% janitor::adorn_pct_formatting(rounding = "half up", 
    digits = 1) %>% knitr::kable()

LVI	n	percent
Absent	148	59.2%
Present	102	40.8%

4.1.2.5 Descriptive Statistics PNI

mydata %>% janitor::tabyl(PNI) %>% janitor::adorn_pct_formatting(rounding = "half up", 
    digits = 1) %>% knitr::kable()

PNI	n	percent	valid_percent
Absent	181	72.4%	72.7%
Present	68	27.2%	27.3%
NA	1	0.4%	-

4.1.2.6 Descriptive Statistics Group

mydata %>% janitor::tabyl(Group) %>% janitor::adorn_pct_formatting(rounding = "half up", 
    digits = 1) %>% knitr::kable()

Group	n	percent	valid_percent
Control	119	47.6%	47.8%
Treatment	130	52.0%	52.2%
NA	1	0.4%	-

4.1.2.7 Descriptive Statistics Grade

mydata %>% janitor::tabyl(Grade) %>% janitor::adorn_pct_formatting(rounding = "half up", 
    digits = 1) %>% knitr::kable()

Grade	n	percent	valid_percent
1	68	27.2%	27.3%
2	84	33.6%	33.7%
3	97	38.8%	39.0%
NA	1	0.4%	-

4.1.2.8 Descriptive Statistics TStage

mydata %>% janitor::tabyl(TStage) %>% janitor::adorn_pct_formatting(rounding = "half up", 
    digits = 1) %>% knitr::kable()

TStage	n	percent
1	23	9.2%
2	49	19.6%
3	65	26.0%
4	113	45.2%

4.1.2.9 Descriptive Statistics LymphNodeMetastasis

mydata %>% janitor::tabyl(LymphNodeMetastasis) %>% janitor::adorn_pct_formatting(rounding = "half up", 
    digits = 1) %>% knitr::kable()

LymphNodeMetastasis	n	percent	valid_percent
Absent	147	58.8%	59.0%
Present	102	40.8%	41.0%
NA	1	0.4%	-

4.1.2.10 Descriptive Statistics Grade_Level

mydata %>% janitor::tabyl(Grade_Level) %>% janitor::adorn_pct_formatting(rounding = "half up", 
    digits = 1) %>% knitr::kable()

Grade_Level	n	percent	valid_percent
high	91	36.4%	36.5%
low	75	30.0%	30.1%
moderate	83	33.2%	33.3%
NA	1	0.4%	-

4.1.2.11 Descriptive Statistics DeathTime

mydata %>% janitor::tabyl(DeathTime) %>% janitor::adorn_pct_formatting(rounding = "half up", 
    digits = 1) %>% knitr::kable()

DeathTime	n	percent
MoreThan1Year	101	40.4%
Within1Year	149	59.6%

race_stats <- summarytools::freq(mydata$Race)
print(race_stats, report.nas = FALSE, totals = FALSE, display.type = FALSE, Variable.label = "Race Group")

mydata %>% explore::describe(PreinvasiveComponent)

variable = PreinvasiveComponent
type     = character
na       = 1 of 250 (0.4%)
unique   = 3
 Absent  = 205 (82%)
 Present = 44 (17.6%)
 NA      = 1 (0.4%)

## Frequency or custom tables for categorical variables
SmartEDA::ExpCTable(mydata, Target = NULL, margin = 1, clim = 10, nlim = 5, round = 2, 
    bin = NULL, per = T)

               Variable         Valid Frequency Percent CumPercent
1                   Sex        Female       122    48.8       48.8
2                   Sex          Male       127    50.8       99.6
3                   Sex            NA         1     0.4      100.0
4                   Sex         TOTAL       250      NA         NA
5                  Race         Asian        11     4.4        4.4
6                  Race     Bi-Racial         5     2.0        6.4
7                  Race         Black        30    12.0       18.4
8                  Race      Hispanic        49    19.6       38.0
9                  Race            NA         1     0.4       38.4
10                 Race        Native         1     0.4       38.8
11                 Race         Other         2     0.8       39.6
12                 Race         White       151    60.4      100.0
13                 Race         TOTAL       250      NA         NA
14 PreinvasiveComponent        Absent       205    82.0       82.0
15 PreinvasiveComponent            NA         1     0.4       82.4
16 PreinvasiveComponent       Present        44    17.6      100.0
17 PreinvasiveComponent         TOTAL       250      NA         NA
18                  LVI        Absent       148    59.2       59.2
19                  LVI       Present       102    40.8      100.0
20                  LVI         TOTAL       250      NA         NA
21                  PNI        Absent       181    72.4       72.4
22                  PNI            NA         1     0.4       72.8
23                  PNI       Present        68    27.2      100.0
24                  PNI         TOTAL       250      NA         NA
25                Group       Control       119    47.6       47.6
26                Group            NA         1     0.4       48.0
27                Group     Treatment       130    52.0      100.0
28                Group         TOTAL       250      NA         NA
29                Grade             1        68    27.2       27.2
30                Grade             2        84    33.6       60.8
31                Grade             3        97    38.8       99.6
32                Grade            NA         1     0.4      100.0
33                Grade         TOTAL       250      NA         NA
34               TStage             1        23     9.2        9.2
35               TStage             2        49    19.6       28.8
36               TStage             3        65    26.0       54.8
37               TStage             4       113    45.2      100.0
38               TStage         TOTAL       250      NA         NA
39  LymphNodeMetastasis        Absent       147    58.8       58.8
40  LymphNodeMetastasis            NA         1     0.4       59.2
41  LymphNodeMetastasis       Present       102    40.8      100.0
42  LymphNodeMetastasis         TOTAL       250      NA         NA
43          Grade_Level          high        91    36.4       36.4
44          Grade_Level           low        75    30.0       66.4
45          Grade_Level      moderate        83    33.2       99.6
46          Grade_Level            NA         1     0.4      100.0
47          Grade_Level         TOTAL       250      NA         NA
48            DeathTime MoreThan1Year       101    40.4       40.4
49            DeathTime   Within1Year       149    59.6      100.0
50            DeathTime         TOTAL       250      NA         NA
51      AntiX_intensity             1        25    10.0       10.0
52      AntiX_intensity             2       103    41.2       51.2
53      AntiX_intensity             3       121    48.4       99.6
54      AntiX_intensity            NA         1     0.4      100.0
55      AntiX_intensity         TOTAL       250      NA         NA
56      AntiY_intensity             1        82    32.8       32.8
57      AntiY_intensity             2       108    43.2       76.0
58      AntiY_intensity             3        59    23.6       99.6
59      AntiY_intensity            NA         1     0.4      100.0
60      AntiY_intensity         TOTAL       250      NA         NA

inspectdf::inspect_cat(mydata)

# A tibble: 16 x 5
   col_name               cnt common      common_pcnt levels            
   <chr>                <int> <chr>             <dbl> <named list>      
 1 Death                    3 TRUE               67.6 <tibble [3 × 3]>  
 2 DeathTime                2 Within1Year        59.6 <tibble [2 × 3]>  
 3 Grade                    4 3                  38.8 <tibble [4 × 3]>  
 4 Grade_Level              4 high               36.4 <tibble [4 × 3]>  
 5 Group                    3 Treatment          52   <tibble [3 × 3]>  
 6 ID                     250 001                 0.4 <tibble [250 × 3]>
 7 LVI                      2 Absent             59.2 <tibble [2 × 3]>  
 8 LymphNodeMetastasis      3 Absent             58.8 <tibble [3 × 3]>  
 9 Name                   250 Aaleyah             0.4 <tibble [250 × 3]>
10 PNI                      3 Absent             72.4 <tibble [3 × 3]>  
11 PreinvasiveComponent     3 Absent             82   <tibble [3 × 3]>  
12 Race                     8 White              60.4 <tibble [8 × 3]>  
13 Sex                      3 Male               50.8 <tibble [3 × 3]>  
14 Smoker                   3 FALSE              52.8 <tibble [3 × 3]>  
15 TStage                   4 4                  45.2 <tibble [4 × 3]>  
16 Valid                    3 FALSE              53.6 <tibble [3 × 3]>  

inspectdf::inspect_cat(mydata)$levels$Group

# A tibble: 3 x 3
  value      prop   cnt
  <chr>     <dbl> <int>
1 Treatment 0.52    130
2 Control   0.476   119
3 <NA>      0.004     1

4.1.2.12 Split-Group Stats Categorical

library(summarytools)

grouped_freqs <- stby(data = mydata$Smoker, INDICES = mydata$Sex, FUN = freq, cumul = FALSE, 
    report.nas = FALSE)

grouped_freqs %>% tb(order = 2)

4.1.2.13 Grouped Categorical

summarytools::stby(list(x = mydata$LVI, y = mydata$LymphNodeMetastasis), mydata$PNI, 
    summarytools::ctable)

with(mydata, summarytools::stby(list(x = LVI, y = LymphNodeMetastasis), PNI, summarytools::ctable))

mydata %>% dplyr::select(characterVariables) %>% dplyr::select(PreinvasiveComponent, 
    PNI, LVI) %>% reactable::reactable(data = ., groupBy = c("PreinvasiveComponent", 
    "PNI"), columns = list(LVI = reactable::colDef(aggregate = "count")))

---

4.1.3 Continious Variables

questionr:::icut()

source(here::here("R", "gc_desc_cont.R"))

Descriptive Statistics Age

mydata %>% jmv::descriptives(data = ., vars = "Age", hist = TRUE, dens = TRUE, box = TRUE, 
    violin = TRUE, dot = TRUE, mode = TRUE, sd = TRUE, variance = TRUE, skew = TRUE, 
    kurt = TRUE, quart = TRUE)

 DESCRIPTIVES

 Descriptives                      
 ───────────────────────────────── 
                          Age      
 ───────────────────────────────── 
   N                         249   
   Missing                     1   
   Mean                     48.6   
   Median                   48.0   
   Mode                     48.0   
   Standard deviation       14.1   
   Variance                  199   
   Minimum                  25.0   
   Maximum                  73.0   
   Skewness               0.0724   
   Std. error skewness     0.154   
   Kurtosis                -1.17   
   Std. error kurtosis     0.307   
   25th percentile          36.0   
   50th percentile          48.0   
   75th percentile          61.0   
 ───────────────────────────────── 

Descriptive Statistics AntiX_intensity

mydata %>% jmv::descriptives(data = ., vars = "AntiX_intensity", hist = TRUE, dens = TRUE, 
    box = TRUE, violin = TRUE, dot = TRUE, mode = TRUE, sd = TRUE, variance = TRUE, 
    skew = TRUE, kurt = TRUE, quart = TRUE)

 DESCRIPTIVES

 Descriptives                               
 ────────────────────────────────────────── 
                          AntiX_intensity   
 ────────────────────────────────────────── 
   N                                  249   
   Missing                              1   
   Mean                              2.39   
   Median                            2.00   
   Mode                              3.00   
   Standard deviation               0.663   
   Variance                         0.439   
   Minimum                           1.00   
   Maximum                           3.00   
   Skewness                        -0.618   
   Std. error skewness              0.154   
   Kurtosis                        -0.648   
   Std. error kurtosis              0.307   
   25th percentile                   2.00   
   50th percentile                   2.00   
   75th percentile                   3.00   
 ────────────────────────────────────────── 

Descriptive Statistics AntiY_intensity

mydata %>% jmv::descriptives(data = ., vars = "AntiY_intensity", hist = TRUE, dens = TRUE, 
    box = TRUE, violin = TRUE, dot = TRUE, mode = TRUE, sd = TRUE, variance = TRUE, 
    skew = TRUE, kurt = TRUE, quart = TRUE)

 DESCRIPTIVES

 Descriptives                               
 ────────────────────────────────────────── 
                          AntiY_intensity   
 ────────────────────────────────────────── 
   N                                  249   
   Missing                              1   
   Mean                              1.91   
   Median                            2.00   
   Mode                              2.00   
   Standard deviation               0.748   
   Variance                         0.560   
   Minimum                           1.00   
   Maximum                           3.00   
   Skewness                         0.152   
   Std. error skewness              0.154   
   Kurtosis                         -1.20   
   Std. error kurtosis              0.307   
   25th percentile                   1.00   
   50th percentile                   2.00   
   75th percentile                   2.00   
 ────────────────────────────────────────── 

tab <- tableone::CreateTableOne(data = mydata)
# ?print.ContTable
tab$ContTable

                             
                              Overall      
  n                           250          
  Age (mean (SD))             48.59 (14.12)
  AntiX_intensity (mean (SD))  2.39 (0.66) 
  AntiY_intensity (mean (SD))  1.91 (0.75) 

print(tab$ContTable, nonnormal = c("Anti-X-intensity"))

                             
                              Overall      
  n                           250          
  Age (mean (SD))             48.59 (14.12)
  AntiX_intensity (mean (SD))  2.39 (0.66) 
  AntiY_intensity (mean (SD))  1.91 (0.75) 

mydata %>% explore::describe(Age)

variable = Age
type     = double
na       = 1 of 250 (0.4%)
unique   = 50
min|max  = 25 | 73
q05|q95  = 27.4 | 71
q25|q75  = 36 | 61
median   = 48
mean     = 48.58635

mydata %>% dplyr::select(continiousVariables) %>% SmartEDA::ExpNumStat(data = ., 
    by = "A", gp = NULL, Qnt = seq(0, 1, 0.1), MesofShape = 2, Outlier = TRUE, round = 2)

inspectdf::inspect_num(mydata, breaks = 10)

# A tibble: 3 x 10
  col_name        min    q1 median  mean    q3   max     sd pcnt_na hist        
  <chr>         <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>  <dbl>   <dbl> <named list>
1 Age              25    36     48 48.6     61    73 14.1       0.4 <tibble [12…
2 AntiX_intens…     1     2      2  2.39     3     3  0.663     0.4 <tibble [12…
3 AntiY_intens…     1     1      2  1.91     2     3  0.748     0.4 <tibble [12…

inspectdf::inspect_num(mydata)$hist$Age

# A tibble: 27 x 2
   value        prop
   <chr>       <dbl>
 1 [-Inf, 24) 0     
 2 [24, 26)   0.0241
 3 [26, 28)   0.0281
 4 [28, 30)   0.0482
 5 [30, 32)   0.0402
 6 [32, 34)   0.0562
 7 [34, 36)   0.0402
 8 [36, 38)   0.0402
 9 [38, 40)   0.0321
10 [40, 42)   0.0402
# … with 17 more rows

inspectdf::inspect_num(mydata, breaks = 10) %>% inspectdf::show_plot()

4.1.3.1 Split-Group Stats Continious

grouped_descr <- summarytools::stby(data = mydata, INDICES = mydata$Sex, FUN = summarytools::descr, 
    stats = "common")
# grouped_descr %>% summarytools::tb(order = 2)
grouped_descr %>% summarytools::tb()

4.1.3.2 Grouped Continious

summarytools::stby(data = mydata, INDICES = mydata$PreinvasiveComponent, FUN = summarytools::descr, 
    stats = c("mean", "sd", "min", "med", "max"), transpose = TRUE)

with(mydata, summarytools::stby(Age, PreinvasiveComponent, summarytools::descr), 
    stats = c("mean", "sd", "min", "med", "max"), transpose = TRUE)

mydata %>% group_by(PreinvasiveComponent) %>% summarytools::descr(stats = "fivenum")

## Summary statistics by – category
SmartEDA::ExpNumStat(mydata, by = "GA", gp = "PreinvasiveComponent", Qnt = seq(0, 
    1, 0.1), MesofShape = 2, Outlier = TRUE, round = 2)

  Vname                        Group  TN nNeg nZero nPos NegInf PosInf NA_Value
1   Age     PreinvasiveComponent:All 250    0     0  249      0      0        1
2   Age PreinvasiveComponent:Present  44    0     0   44      0      0        0
3   Age  PreinvasiveComponent:Absent 205    0     0  204      0      0        1
4   Age      PreinvasiveComponent:NA   0    0     0    0      0      0        0
  Per_of_Missing   sum min  max  mean median    SD   CV   IQR Skewness Kurtosis
1           0.40 12098  25   73 48.59   48.0 14.12 0.29 25.00     0.07    -1.17
2           0.00  2144  25   73 48.73   48.5 15.11 0.31 27.75    -0.05    -1.26
3           0.49  9910  25   73 48.58   48.0 13.96 0.29 25.00     0.10    -1.16
4            NaN     0 Inf -Inf   NaN     NA    NA   NA    NA      NaN      NaN
  0%  10%  20%  30% 40%  50%  60%  70% 80%  90% 100% LB.25% UB.75% nOutliers
1 25 29.8 34.0 39.0  44 48.0 53.0 58.0  64 68.2   73  -1.50  98.50         0
2 25 29.0 30.6 40.0  45 48.5 54.0 58.3  64 69.4   73  -5.38 105.62         0
3 25 31.0 34.0 38.9  44 48.0 52.8 58.0  63 68.0   73  -1.50  98.50         0
4 NA   NA   NA   NA  NA   NA   NA   NA  NA   NA   NA     NA     NA         0

---

4.2 Cross Tables

Codes for cross tables.²⁴

library(finalfit)

dependent <- c("dependent1", "dependent2")

explanatory <- c("explanatory1", "explanatory2")

dependent <- "PreinvasiveComponent"

explanatory <- c("Sex", "Age", "Grade", "TStage")

Change column = TRUE argument to get row or column percentages.

source(here::here("R", "gc_table_cross.R"))

Cross Table PreinvasiveComponent

mydata %>%
    summary_factorlist(dependent = 'PreinvasiveComponent', 
                       explanatory = explanatory,
                       # column = TRUE,
                       total_col = TRUE,
                       p = TRUE,
                       add_dependent_label = TRUE,
                       na_include=FALSE
                       # catTest = catTestfisher
                       ) -> table

knitr::kable(table, row.names = FALSE, align = c('l', 'l', 'r', 'r', 'r'))

Dependent: PreinvasiveComponent		Absent	Present	Total	p
Sex	Female	101 (49.5)	20 (45.5)	121 (48.8)	0.625
	Male	103 (50.5)	24 (54.5)	127 (51.2)
Age	Mean (SD)	48.6 (14.0)	48.7 (15.1)	48.6 (14.1)	0.997
Grade	1	53 (25.9)	15 (34.9)	68 (27.4)	0.453
	2	71 (34.6)	12 (27.9)	83 (33.5)
	3	81 (39.5)	16 (37.2)	97 (39.1)
TStage	1	18 (8.8)	5 (11.4)	23 (9.2)	0.934
	2	41 (20.0)	8 (18.2)	49 (19.7)
	3	52 (25.4)	12 (27.3)	64 (25.7)
	4	94 (45.9)	19 (43.2)	113 (45.4)

library(DT)
datatable(mtcars, rownames = FALSE, filter="top", options = list(pageLength = 5, scrollX=T) )

4.3 Survival Analysis

Codes for Survival Analysis²⁵

• Survival analysis with strata, clusters, frailties and competing risks in in Finalfit

https://www.datasurg.net/2019/09/12/survival-analysis-with-strata-clusters-frailties-and-competing-risks-in-in-finalfit/

• Intracranial WHO grade I meningioma: a competing risk analysis of progression and disease-specific survival

https://link.springer.com/article/10.1007/s00701-019-04096-9

Calculate survival time

mydata$int <- lubridate::interval(lubridate::ymd(mydata$SurgeryDate), lubridate::ymd(mydata$LastFollowUpDate))
mydata$OverallTime <- lubridate::time_length(mydata$int, "month")
mydata$OverallTime <- round(mydata$OverallTime, digits = 1)

recode death status outcome as numbers for survival analysis

## Recoding mydata$Death into mydata$Outcome
mydata$Outcome <- forcats::fct_recode(as.character(mydata$Death), `1` = "TRUE", `0` = "FALSE")
mydata$Outcome <- as.numeric(as.character(mydata$Outcome))

it is always a good practice to double-check after recoding²⁶

table(mydata$Death, mydata$Outcome)

       
          0   1
  FALSE  80   0
  TRUE    0 169

4.3.1 Kaplan-Meier

library(survival)
# data(lung) km <- with(lung, Surv(time, status))
km <- with(mydata, Surv(OverallTime, Outcome))
head(km, 80)

 [1] 10.5   9.1   7.1+  8.0   5.8   8.2+  8.8+  9.3   9.0  10.5+  3.4+  8.7 
[13]  3.8   5.6   8.1   8.9   8.1   8.9   7.7+  5.2+ 10.5   4.0   7.2   5.2+
[25]  4.2   3.3   4.0   4.0+ 11.0  10.7+  7.3+ 11.4+  6.4   8.7    NA+  6.9+
[37]  5.9  10.5+  8.8   8.2  11.5+  7.0   6.4   9.6   5.6+  9.9+  8.2+  9.1 
[49]  3.3  11.6   7.3+ 10.8   7.6  11.5   4.0+  8.1+ 10.4+  6.2+  9.6   6.5 
[61] 11.8+  4.1   6.6   5.8+  9.3+  5.4   7.1   9.6   3.9   4.9   3.3   9.0 
[73]  5.8   9.3   5.2   8.3  10.3  10.6   7.1+  6.4 

plot(km)

Kaplan-Meier Plot Log-Rank Test

# Drawing Survival Curves Using ggplot2
# https://rpkgs.datanovia.com/survminer/reference/ggsurvplot.html
dependentKM <- "Surv(OverallTime, Outcome)"
explanatoryKM <- "LVI"

mydata %>%
  finalfit::surv_plot(.data = .,
                      dependent = dependentKM,
                      explanatory = explanatoryKM,
                      xlab='Time (months)',
                      pval=TRUE,
                      legend = 'none',
                      break.time.by = 12,
                      xlim = c(0,60)
                      # legend.labs = c('a','b')
                      )

# Drawing Survival Curves Using ggplot2
# https://rpkgs.datanovia.com/survminer/reference/ggsurvplot.html

mydata %>%
  finalfit::surv_plot(.data = .,
                      dependent = "Surv(OverallTime, Outcome)",
                      explanatory = "LVI",
                      xlab='Time (months)',
                      pval=TRUE,
                      legend = 'none',
                      break.time.by = 12,
                      xlim = c(0,60)
                      # legend.labs = c('a','b')
                      )

4.3.2 Univariate Cox-Regression

library(finalfit)
library(survival)
explanatoryUni <- "LVI"
dependentUni <- "Surv(OverallTime, Outcome)"

tUni <- mydata %>% finalfit::finalfit(dependentUni, explanatoryUni)

knitr::kable(tUni, row.names = FALSE, align = c("l", "l", "r", "r", "r", "r"))

Dependent: Surv(OverallTime, Outcome)		all	HR (univariable)	HR (multivariable)
LVI	Absent	148 (100.0)	-	-
	Present	102 (100.0)	1.45 (1.06-1.98, p=0.021)	1.45 (1.06-1.98, p=0.021)

tUni_df <- tibble::as_tibble(tUni, .name_repair = "minimal") %>% janitor::clean_names()

tUni_df_descr <- paste0("When ", tUni_df$dependent_surv_overall_time_outcome[1], 
    " is ", tUni_df$x[2], ", there is ", tUni_df$hr_univariable[2], " times risk than ", 
    "when ", tUni_df$dependent_surv_overall_time_outcome[1], " is ", tUni_df$x[1], 
    ".")

When LVI is Present, there is 1.45 (1.06-1.98, p=0.021) times risk than when LVI is Absent.

4.3.3 Kaplan-Meier Median Survival

km_fit <- survfit(Surv(OverallTime, Outcome) ~ LVI, data = mydata)
km_fit

Call: survfit(formula = Surv(OverallTime, Outcome) ~ LVI, data = mydata)

   3 observations deleted due to missingness 
              n events median 0.95LCL 0.95UCL
LVI=Absent  146    100   19.7    14.1    27.1
LVI=Present 101     69   10.5     9.3    12.8

plot(km_fit)

# summary(km_fit)

km_fit_median_df <- summary(km_fit)
km_fit_median_df <- as.data.frame(km_fit_median_df$table) %>% janitor::clean_names() %>% 
    tibble::rownames_to_column()

km_fit_median_definition <- km_fit_median_df %>% dplyr::mutate(description = glue::glue("When {rowname}, median survival is {median} [{x0_95lcl} - {x0_95ucl}, 95% CI] months.")) %>% 
    dplyr::select(description) %>% dplyr::pull()

When LVI=Absent, median survival is 19.7 [14.1 - 27.1, 95% CI] months., When LVI=Present, median survival is 10.5 [9.3 - 12.8, 95% CI] months.

4.3.4 1-3-5-yr survival

summary(km_fit, times = c(12, 36, 60))

Call: survfit(formula = Surv(OverallTime, Outcome) ~ LVI, data = mydata)

3 observations deleted due to missingness 
                LVI=Absent 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12     71      52    0.612  0.0426        0.534        0.701
   36     19      38    0.235  0.0424        0.165        0.334

                LVI=Present 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12     27      52    0.391  0.0544       0.2978        0.514
   36      5      15    0.120  0.0436       0.0593        0.245

km_fit_summary <- summary(km_fit, times = c(12, 36, 60))

km_fit_df <- as.data.frame(km_fit_summary[c("strata", "time", "n.risk", "n.event", 
    "surv", "std.err", "lower", "upper")])

km_fit_definition <- km_fit_df %>% dplyr::mutate(description = glue::glue("When {strata}, {time} month survival is {scales::percent(surv)} [{scales::percent(lower)}-{scales::percent(upper)}, 95% CI].")) %>% 
    dplyr::select(description) %>% dplyr::pull()

When LVI=Absent, 12 month survival is 61% [53%-70.1%, 95% CI]., When LVI=Absent, 36 month survival is 23% [16%-33.4%, 95% CI]., When LVI=Present, 12 month survival is 39% [30%-51.4%, 95% CI]., When LVI=Present, 36 month survival is 12% [6%-24.5%, 95% CI].

source(here::here("R", "gc_survival.R"))

4.3.5 Survival Analysis LVI

Kaplan-Meier Plot Log-Rank Test

library(survival)
library(survminer)
library(finalfit)

mydata %>%
  finalfit::surv_plot('Surv(OverallTime, Outcome)', 'LVI', 
  xlab='Time (months)', pval=TRUE, legend = 'none',
    break.time.by = 12, xlim = c(0,60)

# legend.labs = c('a','b')

)

Univariate Cox-Regression

explanatoryUni <- "LVI"
dependentUni <- "Surv(OverallTime, Outcome)"
tUni <- mydata %>% finalfit(dependentUni, explanatoryUni)

knitr::kable(tUni, row.names = FALSE, align = c("l", "l", "r", "r", "r", "r"))

Dependent: Surv(OverallTime, Outcome)		all	HR (univariable)	HR (multivariable)
LVI	Absent	148 (100.0)	-	-
	Present	102 (100.0)	1.45 (1.06-1.98, p=0.021)	1.45 (1.06-1.98, p=0.021)

Univariate Cox-Regression Summary

tUni_df <- tibble::as_tibble(tUni, .name_repair = "minimal") %>% janitor::clean_names(dat = ., 
    case = "snake")


n_level <- dim(tUni_df)[1]

tUni_df_descr <- function(n) {
    paste0("When ", tUni_df$dependent_surv_overall_time_outcome[1], " is ", tUni_df$x[n + 
        1], ", there is ", tUni_df$hr_univariable[n + 1], " times risk than ", "when ", 
        tUni_df$dependent_surv_overall_time_outcome[1], " is ", tUni_df$x[1], ".")
    
}



results5 <- purrr::map(.x = c(2:n_level - 1), .f = tUni_df_descr)

print(unlist(results5))

[1] "When LVI is Present, there is 1.45 (1.06-1.98, p=0.021) times risk than when LVI is Absent."

Median Survival

km_fit <- survfit(Surv(OverallTime, Outcome) ~ LVI, data = mydata)
km_fit

Call: survfit(formula = Surv(OverallTime, Outcome) ~ LVI, data = mydata)

   3 observations deleted due to missingness 
              n events median 0.95LCL 0.95UCL
LVI=Absent  146    100   19.7    14.1    27.1
LVI=Present 101     69   10.5     9.3    12.8

km_fit_median_df <- summary(km_fit)
km_fit_median_df <- as.data.frame(km_fit_median_df$table) %>% janitor::clean_names(dat = ., 
    case = "snake") %>% tibble::rownames_to_column(.data = ., var = "LVI")



km_fit_median_definition <- km_fit_median_df %>% dplyr::mutate(description = glue::glue("When, LVI, {LVI}, median survival is {median} [{x0_95lcl} - {x0_95ucl}, 95% CI] months.")) %>% 
    dplyr::mutate(description = gsub(pattern = "thefactor=", replacement = " is ", 
        x = description)) %>% dplyr::select(description) %>% dplyr::pull()

km_fit_median_definition

When, LVI, LVI=Absent, median survival is 19.7 [14.1 - 27.1, 95% CI] months.
When, LVI, LVI=Present, median survival is 10.5 [9.3 - 12.8, 95% CI] months.

1-3-5-yr survival

summary(km_fit, times = c(12, 36, 60))

Call: survfit(formula = Surv(OverallTime, Outcome) ~ LVI, data = mydata)

3 observations deleted due to missingness 
                LVI=Absent 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12     71      52    0.612  0.0426        0.534        0.701
   36     19      38    0.235  0.0424        0.165        0.334

                LVI=Present 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12     27      52    0.391  0.0544       0.2978        0.514
   36      5      15    0.120  0.0436       0.0593        0.245

km_fit_summary <- summary(km_fit, times = c(12, 36, 60))

km_fit_df <- as.data.frame(km_fit_summary[c("strata", "time", "n.risk", "n.event", 
    "surv", "std.err", "lower", "upper")])

km_fit_df

       strata time n.risk n.event      surv    std.err      lower     upper
1  LVI=Absent   12     71      52 0.6115055 0.04257648 0.53350053 0.7009157
2  LVI=Absent   36     19      38 0.2347686 0.04235846 0.16483932 0.3343638
3 LVI=Present   12     27      52 0.3911397 0.05442593 0.29777603 0.5137763
4 LVI=Present   36      5      15 0.1204428 0.04355881 0.05928425 0.2446936

km_fit_definition <- km_fit_df %>% dplyr::mutate(description = glue::glue("When {strata}, {time} month survival is {scales::percent(surv)} [{scales::percent(lower)}-{scales::percent(upper)}, 95% CI].")) %>% 
    dplyr::select(description) %>% dplyr::pull()

km_fit_definition

When LVI=Absent, 12 month survival is 61% [53%-70.1%, 95% CI].
When LVI=Absent, 36 month survival is 23% [16%-33.4%, 95% CI].
When LVI=Present, 12 month survival is 39% [30%-51.4%, 95% CI].
When LVI=Present, 36 month survival is 12% [6%-24.5%, 95% CI].

summary(km_fit)$table

            records n.max n.start events   *rmean *se(rmean) median 0.95LCL
LVI=Absent      146   146     146    100 24.32960   1.653307   19.7    14.1
LVI=Present     101   101     101     69 18.16069   1.840202   10.5     9.3
            0.95UCL
LVI=Absent     27.1
LVI=Present    12.8

km_fit_median_df <- summary(km_fit)
results1html <- as.data.frame(km_fit_median_df$table) %>% janitor::clean_names(dat = ., 
    case = "snake") %>% tibble::rownames_to_column(.data = ., var = "LVI")

results1html[, 1] <- gsub(pattern = "thefactor=", replacement = "", x = results1html[, 
    1])

knitr::kable(results1html, row.names = FALSE, align = c("l", rep("r", 9)), format = "html", 
    digits = 1)

LVI	records	n_max	n_start	events	rmean	se_rmean	median	x0_95lcl	x0_95ucl
LVI=Absent	146	146	146	100	24.3	1.7	19.7	14.1	27.1
LVI=Present	101	101	101	69	18.2	1.8	10.5	9.3	12.8

Pairwise Comparisons

4.3.6 Pairwise comparison

dependentKM <- "Surv(OverallTime, Outcome)"
explanatoryKM <- "TStage"

mydata %>%
  finalfit::surv_plot(.data = .,
                      dependent = dependentKM,
                      explanatory = explanatoryKM,
                      xlab='Time (months)',
                      pval=TRUE,
                      legend = 'none',
                      break.time.by = 12,
                      xlim = c(0,60)
                      # legend.labs = c('a','b')
                      )

km_fit

Call: survfit(formula = Surv(OverallTime, Outcome) ~ LVI, data = mydata)

   3 observations deleted due to missingness 
              n events median 0.95LCL 0.95UCL
LVI=Absent  146    100   19.7    14.1    27.1
LVI=Present 101     69   10.5     9.3    12.8

print(km_fit, 
      scale=1,
      digits = max(options()$digits - 4,3),
      print.rmean=getOption("survfit.print.rmean"),
      rmean = getOption('survfit.rmean'),
      print.median=getOption("survfit.print.median"),
      median = getOption('survfit.median')
      
      )

Call: survfit(formula = Surv(OverallTime, Outcome) ~ LVI, data = mydata)

   3 observations deleted due to missingness 
              n events median 0.95LCL 0.95UCL
LVI=Absent  146    100   19.7    14.1    27.1
LVI=Present 101     69   10.5     9.3    12.8

4.3.7 Multivariate Analysis Survival

---

---

5 Discussion

• Interpret the results in context of the working hypothesis elaborated in the introduction and other relevant studies; include a discussion of limitations of the study.

• Discuss potential clinical applications and implications for future research

6 Footer

Codes for explaining the software and the packages that are used in the analysis²⁷

6.0.1 Save Final Data

projectName <- list.files(path = here::here(), pattern = "Rproj")
projectName <- gsub(pattern = ".Rproj", replacement = "", x = projectName)

analysisDate <- as.character(Sys.Date())

imageName <- paste0(projectName, analysisDate, ".RData")

save.image(file = here::here("data", imageName))

rdsName <- paste0(projectName, analysisDate, ".rds")

readr::write_rds(x = mydata, path = here::here("data", rdsName))

saveRDS(object = mydata, file = here::here("data", rdsName))

excelName <- paste0(projectName, analysisDate, ".xlsx")

rio::export(x = mydata, file = here::here("data", excelName), format = "xlsx")

# writexl::write_xlsx(mydata, here::here('data', excelName))

print(glue::glue("saved data after analysis to ", rownames(file.info(here::here("data", 
    excelName))), " : ", as.character(file.info(here::here("data", excelName))$ctime)))

saved data after analysis to /Users/serdarbalciold/histopathRprojects/histopathology-template/data/histopathology-template2020-02-18.xlsx : 2020-02-18 18:47:50

mydata %>% downloadthis::download_this(output_name = excelName, output_extension = ".csv", 
    button_label = "Download data as csv", button_type = "default")

Download data as csv

mydata %>% downloadthis::download_this(output_name = excelName, output_extension = ".xlsx", 
    button_label = "Download data as xlsx", button_type = "primary")

Download data as xlsx

---

6.0.2 Final Data Summary

# use summarytools to generate final data summary
# summarytools::view(summarytools::dfSummary(x = mydata
#                                            , style = "markdown"))

---

6.0.3 Software and Libraries Used

Why and how to cite software and packages?²⁸

citation()

To cite R in publications use:

  R Core Team (2019). R: A language and environment for statistical
  computing. R Foundation for Statistical Computing, Vienna, Austria.
  URL https://www.R-project.org/.

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {R: A Language and Environment for Statistical Computing},
    author = {{R Core Team}},
    organization = {R Foundation for Statistical Computing},
    address = {Vienna, Austria},
    year = {2019},
    url = {https://www.R-project.org/},
  }

We have invested a lot of time and effort in creating R, please cite it
when using it for data analysis. See also 'citation("pkgname")' for
citing R packages.

The jamovi project (2019). jamovi. (Version 0.9) [Computer Software]. Retrieved from https://www.jamovi.org.
R Core Team (2018). R: A Language and envionment for statistical computing. [Computer software]. Retrieved from https://cran.r-project.org/.
Fox, J., & Weisberg, S. (2018). car: Companion to Applied Regression. [R package]. Retrieved from https://cran.r-project.org/package=car. Wickham et al., (2019). Welcome to the tidyverse. Journal of Open Source Software, 4(43), 1686, https://doi.org/10.21105/joss.01686 Data processing was carried out with R (R Core Team, 2019) and the easystats ecosystem (Lüdecke, Waggoner, & Makowski, 2019; Makowski, Ben-Shachar, & Lüdecke, 2019)

report::cite_packages(session = sessionInfo())

6.1 References

Alastair Rushworth (2019). inspectdf: Inspection, Comparison and Visualisation of Data Frames. R package version 0.0.7. https://CRAN.R-project.org/package=inspectdf
Alboukadel Kassambara (2020). ggpubr: ‘ggplot2’ Based Publication Ready Plots. R package version 0.2.5. https://CRAN.R-project.org/package=ggpubr
Alboukadel Kassambara, Marcin Kosinski and Przemyslaw Biecek (2019). survminer: Drawing Survival Curves using ‘ggplot2’. R package version 0.4.6. https://CRAN.R-project.org/package=survminer
Benjamin Elbers (2020). tidylog: Logging for ‘dplyr’ and ‘tidyr’ Functions. R package version 1.0.0. https://CRAN.R-project.org/package=tidylog
Boxuan Cui (2020). DataExplorer: Automate Data Exploration and Treatment. R package version 0.8.1. https://CRAN.R-project.org/package=DataExplorer
Chung-hong Chan, Geoffrey CH Chan, Thomas J. Leeper, and Jason Becker (2018). rio: A Swiss-army knife for data file I/O. R package version 0.5.16.
David Robinson and Alex Hayes (2020). broom: Convert Statistical Analysis Objects into Tidy Tibbles. R package version 0.5.4. https://CRAN.R-project.org/package=broom
Dayanand Ubrangala, Kiran R, Ravi Prasad Kondapalli and Sayan Putatunda (2020). SmartEDA: Summarize and Explore the Data. R package version 0.3.3. https://CRAN.R-project.org/package=SmartEDA
Dirk Eddelbuettel and Romain Francois (2011). Rcpp: Seamless R and C++ Integration. Journal of Statistical Software, 40(8), 1-18. URL http://www.jstatsoft.org/v40/i08/.
Dirk Eddelbuettel with contributions by Antoine Lucas, Jarek Tuszynski, Henrik Bengtsson, Simon Urbanek, Mario Frasca, Bryan Lewis, Murray Stokely, Hannes Muehleisen, Duncan Murdoch, Jim Hester, Wush Wu, Qiang Kou, Thierry Onkelinx, Michel Lang, Viliam Simko, Kurt Hornik, Radford Neal, Kendon Bell, Matthew de Queljoe, Ion Suruceanu and Bill Denney. (2020). digest: Create Compact Hash Digests of R Objects. R package version 0.6.24. https://CRAN.R-project.org/package=digest Ethan Heinzen, Jason Sinnwell, Elizabeth Atkinson, Tina Gunderson and Gregory Dougherty (2020). arsenal: An Arsenal of ‘R’ Functions for Large-Scale Statistical Summaries. R package version 3.4.0. https://CRAN.R-project.org/package=arsenal
Ewen Harrison, Tom Drake and Riinu Ots (2019). finalfit: Quickly Create Elegant Regression Results Tables and Plots when Modelling. R package version 0.9.7. https://CRAN.R-project.org/package=finalfit
Garrett Grolemund, Hadley Wickham (2011). Dates and Times Made Easy with lubridate. Journal of Statistical Software, 40(3), 1-25. URL http://www.jstatsoft.org/v40/i03/.
H. Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016.
Hadley Wickham (2019). feather: R Bindings to the Feather ‘API’. R package version 0.3.5. https://CRAN.R-project.org/package=feather
Hadley Wickham (2019). forcats: Tools for Working with Categorical Variables (Factors). R package version 0.4.0. https://CRAN.R-project.org/package=forcats
Hadley Wickham (2019). httr: Tools for Working with URLs and HTTP. R package version 1.4.1. https://CRAN.R-project.org/package=httr
Hadley Wickham (2019). modelr: Modelling Functions that Work with the Pipe. R package version 0.1.5. https://CRAN.R-project.org/package=modelr
Hadley Wickham (2019). rvest: Easily Harvest (Scrape) Web Pages. R package version 0.3.5. https://CRAN.R-project.org/package=rvest
Hadley Wickham (2019). stringr: Simple, Consistent Wrappers for Common String Operations. R package version 1.4.0. https://CRAN.R-project.org/package=stringr
Hadley Wickham and Evan Miller (2019). haven: Import and Export ‘SPSS’, ‘Stata’ and ‘SAS’ Files. R package version 2.2.0. https://CRAN.R-project.org/package=haven
Hadley Wickham and Jennifer Bryan (2019). readxl: Read Excel Files. R package version 1.3.1. https://CRAN.R-project.org/package=readxl
Hadley Wickham and Lionel Henry (2020). tidyr: Tidy Messy Data. R package version 1.0.2. https://CRAN.R-project.org/package=tidyr
Hadley Wickham and Yihui Xie (2019). evaluate: Parsing and Evaluation Tools that Provide More Details than the Default. R package version 0.14. https://CRAN.R-project.org/package=evaluate
Hadley Wickham, Jim Hester and Jeroen Ooms (2019). xml2: Parse XML. R package version 1.2.2. https://CRAN.R-project.org/package=xml2
Hadley Wickham, Jim Hester and Romain Francois (2018). readr: Read Rectangular Text Data. R package version 1.3.1. https://CRAN.R-project.org/package=readr
Hadley Wickham, Romain François, Lionel Henry and Kirill Müller (2020). dplyr: A Grammar of Data Manipulation. R package version 0.8.4. https://CRAN.R-project.org/package=dplyr
Jeremy Stephens, Kirill Simonov, Yihui Xie, Zhuoer Dong, Hadley Wickham, Jeffrey Horner, reikoch, Will Beasley, Brendan O’Connor and Gregory R. Warnes (2020). yaml: Methods to Convert R Data to YAML and Back. R package version 2.2.1. https://CRAN.R-project.org/package=yaml
Jeroen Ooms (2014). The jsonlite Package: A Practical and Consistent Mapping Between JSON Data and R Objects. arXiv:1403.2805 [stat.CO] URL https://arxiv.org/abs/1403.2805.
Jim Hester (2019). glue: Interpreted String Literals. R package version 1.3.1. https://CRAN.R-project.org/package=glue
Jim Hester and Gábor Csárdi (2019). pak: Another Approach to Package Installation. R package version 0.1.2. https://CRAN.R-project.org/package=pak
Jim Hester, Gábor Csárdi, Hadley Wickham, Winston Chang, Martin Morgan and Dan Tenenbaum (2020). remotes: R Package Installation from Remote Repositories, Including ‘GitHub’. R package version 2.1.1. https://CRAN.R-project.org/package=remotes
JJ Allaire and Yihui Xie and Jonathan McPherson and Javier Luraschi and Kevin Ushey and Aron Atkins and Hadley Wickham and Joe Cheng and Winston Chang and Richard Iannone (2020). rmarkdown: Dynamic Documents for R. R package version 2.1. URL https://rmarkdown.rstudio.com.
JJ Allaire, Jeffrey Horner, Yihui Xie, Vicent Marti and Natacha Porte (2019). markdown: Render Markdown with the C Library ‘Sundown’. R package version 1.1. https://CRAN.R-project.org/package=markdown
Kazuki Yoshida (2019). tableone: Create ‘Table 1’ to Describe Baseline Characteristics. R package version 0.10.0. https://CRAN.R-project.org/package=tableone
Kevin Ushey (2020). renv: Project Environments. R package version 0.9.3. https://CRAN.R-project.org/package=renv
Kirill Müller (2017). here: A Simpler Way to Find Your Files. R package version 0.1. https://CRAN.R-project.org/package=here
Kirill Müller (2020). hms: Pretty Time of Day. R package version 0.5.3. https://CRAN.R-project.org/package=hms
Kirill Müller and Hadley Wickham (2019). tibble: Simple Data Frames. R package version 2.1.3. https://CRAN.R-project.org/package=tibble
Koji Makiyama (2016). magicfor: Magic Functions to Obtain Results from for Loops. R package version 0.1.0. https://CRAN.R-project.org/package=magicfor
Lionel Henry and Hadley Wickham (2019). purrr: Functional Programming Tools. R package version 0.3.3. https://CRAN.R-project.org/package=purrr
Lionel Henry and Hadley Wickham (2020). rlang: Functions for Base Types and Core R and ‘Tidyverse’ Features. R package version 0.4.4. https://CRAN.R-project.org/package=rlang
Makowski, D. & Lüdecke, D. (2019). The report package for R: Ensuring the use of best practices for results reporting. CRAN. Available from https://github.com/easystats/report. doi: .
Pablo Seibelt (2017). xray: X Ray Vision on your Datasets. R package version 0.2. https://CRAN.R-project.org/package=xray
Paul Hendricks (2015). describer: Describe Data in R Using Common Descriptive Statistics. R package version 0.2.0. https://CRAN.R-project.org/package=describer
Petersen AH, Ekstrøm CT (2019). “dataMaid: Your Assistant forDocumenting Supervised Data Quality Screening in R.” Journal ofStatistical Software, 90(6), 1-38. doi: 10.18637/jss.v090.i06 (URL:https://doi.org/10.18637/jss.v090.i06).
Rinker, T. W. (2018). wakefield: Generate Random Data. version 0.3.3. Buffalo, New York. https://github.com/trinker/wakefield
Rinker, T. W. & Kurkiewicz, D. (2017). pacman: Package Management for R. version 0.5.0. Buffalo, New York. http://github.com/trinker/pacman
Roland Krasser (2020). explore: Simplifies Exploratory Data Analysis. R package version 0.5.4. https://CRAN.R-project.org/package=explore
RStudio and Inc. (2019). htmltools: Tools for HTML. R package version 0.4.0. https://CRAN.R-project.org/package=htmltools
Sam Firke (2020). janitor: Simple Tools for Examining and Cleaning Dirty Data. R package version 1.2.1. https://CRAN.R-project.org/package=janitor
Simon Garnier (2018). viridis: Default Color Maps from ‘matplotlib’. R package version 0.5.1. https://CRAN.R-project.org/package=viridis
Simon Garnier (2018). viridisLite: Default Color Maps from ‘matplotlib’ (Lite Version). R package version 0.3.0. https://CRAN.R-project.org/package=viridisLite
Simon Urbanek (2015). base64enc: Tools for base64 encoding. R package version 0.1-3. https://CRAN.R-project.org/package=base64enc
Stefan Milton Bache and Hadley Wickham (2014). magrittr: A Forward-Pipe Operator for R. R package version 1.5. https://CRAN.R-project.org/package=magrittr
Therneau T (2015). A Package for Survival Analysis in S. version2.38, <URL: https://CRAN.R-project.org/package=survival>.
Tierney N (2017). “visdat: Visualising Whole Data Frames.” JOSS,2(16), 355. doi: 10.21105/joss.00355 (URL:https://doi.org/10.21105/joss.00355), <URL:http://dx.doi.org/10.21105/joss.00355>.
Winston Chang, Joe Cheng, JJ Allaire, Yihui Xie and Jonathan McPherson (2019). shiny: Web Application Framework for R. R package version 1.4.0. https://CRAN.R-project.org/package=shiny
Yihui Xie (2019). formatR: Format R Code Automatically. R package version 1.7. https://CRAN.R-project.org/package=formatR
Yihui Xie (2020). knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.28.
Yihui Xie (2020). mime: Map Filenames to MIME Types. R package version 0.9. https://CRAN.R-project.org/package=mime
Yixuan Qiu and Yihui Xie (2019). highr: Syntax Highlighting for R Source Code. R package version 0.8. https://CRAN.R-project.org/package=highr

report::show_packages(session = sessionInfo()) %>% kableExtra::kable()

# citation('tidyverse')
citation("readxl")

To cite package 'readxl' in publications use:

  Hadley Wickham and Jennifer Bryan (2019). readxl: Read Excel Files. R
  package version 1.3.1. https://CRAN.R-project.org/package=readxl

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {readxl: Read Excel Files},
    author = {Hadley Wickham and Jennifer Bryan},
    year = {2019},
    note = {R package version 1.3.1},
    url = {https://CRAN.R-project.org/package=readxl},
  }

citation("janitor")

To cite package 'janitor' in publications use:

  Sam Firke (2020). janitor: Simple Tools for Examining and Cleaning
  Dirty Data. R package version 1.2.1.
  https://CRAN.R-project.org/package=janitor

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {janitor: Simple Tools for Examining and Cleaning Dirty Data},
    author = {Sam Firke},
    year = {2020},
    note = {R package version 1.2.1},
    url = {https://CRAN.R-project.org/package=janitor},
  }

# citation('report')
citation("finalfit")

To cite package 'finalfit' in publications use:

  Ewen Harrison, Tom Drake and Riinu Ots (2019). finalfit: Quickly
  Create Elegant Regression Results Tables and Plots when Modelling. R
  package version 0.9.7. https://CRAN.R-project.org/package=finalfit

A BibTeX entry for LaTeX users is

  @Manual{,
    title = {finalfit: Quickly Create Elegant Regression Results Tables and Plots when
Modelling},
    author = {Ewen Harrison and Tom Drake and Riinu Ots},
    year = {2019},
    note = {R package version 0.9.7},
    url = {https://CRAN.R-project.org/package=finalfit},
  }

# citation('ggstatsplot')

if (!dir.exists(here::here("bib"))) {
    dir.create(here::here("bib"))
}

knitr::write_bib(x = c(.packages(), "knitr", "shiny"), file = here::here("bib", "packages.bib"))

---

6.1.1 Session Info

sessionInfo()

R version 3.6.0 (2019-04-26)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.15.3

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices datasets  utils     methods   base     

other attached packages:
 [1] survminer_0.4.6    ggpubr_0.2.5       viridis_0.5.1      viridisLite_0.3.0 
 [5] shiny_1.4.0        survival_3.1-8     magrittr_1.5       report_0.1.0      
 [9] wakefield_0.3.3    SmartEDA_0.3.3     magicfor_0.1.0     tableone_0.10.0   
[13] arsenal_3.4.0      DataExplorer_0.8.1 xray_0.2           visdat_0.5.3      
[17] inspectdf_0.0.7    describer_0.2.0    dataMaid_1.4.0     finalfit_0.9.7    
[21] explore_0.5.4      rio_0.5.16         janitor_1.2.1      formatR_1.7       
[25] renv_0.9.3         rlang_0.4.4        glue_1.3.1         tidylog_1.0.0     
[29] broom_0.5.4        modelr_0.1.5       rvest_0.3.5        xml2_1.2.2        
[33] readxl_1.3.1       httr_1.4.1         haven_2.2.0        feather_0.3.5     
[37] lubridate_1.7.4    hms_0.5.3          forcats_0.4.0      stringr_1.4.0     
[41] tibble_2.1.3       purrr_0.3.3        readr_1.3.1        tidyr_1.0.2       
[45] dplyr_0.8.4        ggplot2_3.2.1      rmarkdown_2.1      mime_0.9          
[49] base64enc_0.1-3    jsonlite_1.6.1     knitr_1.28         htmltools_0.4.0   
[53] Rcpp_1.0.3         yaml_2.2.1         markdown_1.1       highr_0.8         
[57] digest_0.6.24      evaluate_0.14      here_0.1           pak_0.1.2         
[61] pacman_0.5.1       remotes_2.1.1     

loaded via a namespace (and not attached):
  [1] utf8_1.1.4          tidyselect_1.0.0    lme4_1.1-21        
  [4] htmlwidgets_1.5.1   grid_3.6.0          lpSolve_5.6.15     
  [7] munsell_0.5.0       effectsize_0.1.2    codetools_0.2-16   
 [10] DT_0.12.1           withr_2.1.2         ISLR_1.2           
 [13] colorspace_1.4-1    rstudioapi_0.11     robustbase_0.93-5  
 [16] ggsignif_0.6.0      labeling_0.3        KMsurv_0.1-5       
 [19] farver_2.0.3        rprojroot_1.3-2     vctrs_0.2.2        
 [22] generics_0.0.2      xfun_0.12           downloadthis_0.1.0 
 [25] R6_2.4.1            reshape_0.8.8       assertthat_0.2.1   
 [28] promises_1.1.0      networkD3_0.4       scales_1.1.0       
 [31] nnet_7.3-12         gtable_0.3.0        clisymbols_1.2.0   
 [34] whoami_1.3.0        splines_3.6.0       lazyeval_0.2.2     
 [37] acepack_1.4.1       bsplus_0.1.1        checkmate_2.0.0    
 [40] backports_1.1.5     httpuv_1.5.2        Hmisc_4.3-1        
 [43] tools_3.6.0         ellipsis_0.3.0      RColorBrewer_1.1-2 
 [46] plyr_1.8.5          jmvcore_1.2.5       progress_1.2.2     
 [49] prettyunits_1.1.1   rpart_4.1-15        sampling_2.8       
 [52] zoo_1.8-7           reactR_0.4.2        cluster_2.1.0      
 [55] fs_1.3.1            survey_3.37         data.table_1.12.8  
 [58] openxlsx_4.1.4      mitml_0.3-7         reactable_0.1.0    
 [61] xtable_1.8-4        jpeg_0.1-8.1        gridExtra_2.3      
 [64] compiler_3.6.0      mice_3.7.0          writexl_1.2        
 [67] crayon_1.3.4        minqa_1.2.4         later_1.0.0        
 [70] Formula_1.2-3       DBI_1.1.0           jmv_1.2.5          
 [73] MASS_7.3-51.5       boot_1.3-24         Matrix_1.2-18      
 [76] cli_2.0.1           mitools_2.4         parallel_3.6.0     
 [79] insight_0.8.1       pan_1.6             igraph_1.2.4.2     
 [82] pkgconfig_2.0.3     km.ci_0.5-2         foreign_0.8-75     
 [85] foreach_1.4.8       snakecase_0.11.0    parameters_0.5.0.1 
 [88] cellranger_1.1.0    survMisc_0.5.5      htmlTable_1.13.3   
 [91] curl_4.3            jomo_2.6-10         rjson_0.2.20       
 [94] nloptr_1.2.1        lifecycle_0.1.0     nlme_3.1-144       
 [97] fansi_0.4.1         labelled_2.2.2      pillar_1.4.3       
[100] ggsci_2.9           lattice_0.20-38     GGally_1.4.0       
[103] fastmap_1.0.1       DEoptimR_1.0-8      bayestestR_0.5.2   
[106] zip_2.0.4           png_0.1-7           iterators_1.0.12   
[109] pander_0.6.3        performance_0.4.4   class_7.3-15       
[112] stringi_1.4.6       ggfittext_0.8.1     latticeExtra_0.6-29
[115] e1071_1.7-3        

---

6.1.2 Loaded packages

pacman::p_loaded(all = TRUE)

---

6.1.3 Notes

Last update on 2020-02-18 18:47:55

Serdar Balci, MD, Pathologist
drserdarbalci@gmail.com
https://rpubs.com/sbalci/CV

---

6.2 Code Appendix

Use following chunk options to include all codes below the report.

{r, echo=TRUE, eval=FALSE, ref.label=knitr::all_labels()}

# installing necessary packages
if (requireNamespace("magrittr", quietly = TRUE)) {
  `%>%` <- magrittr::`%>%`
}
if (!require("remotes")) install.packages("remotes")
if (!require("pacman")) install.packages("pacman")
if (!require("pak")) install.packages("pak")
if (!require("here")) install.packages("here")
source_rmd <- function(rmd_file){
  knitr::knit(rmd_file, output = tempfile(), envir = globalenv())
}

list_of_Rmd <- list.files(path = here::here("childRmd"), pattern = "Rmd")

list_of_Rmd <- list_of_Rmd[!list_of_Rmd %in% c("_19shinySurvival.Rmd")]

purrr::map(.x = here::here("childRmd", list_of_Rmd), .f = source_rmd)

knitr::opts_chunk$set(
    eval = TRUE,
    echo = TRUE,
    fig.path = here::here("figs/"),
    message = FALSE,
    warning = FALSE,
    error = FALSE,
    cache = TRUE,
    comment = NA,
    tidy = TRUE,
    fig.width = 6,
    fig.height = 4
)
# https://cran.r-project.org/web/packages/exploreR/vignettes/exploreR.html
# exploreR::reset()
source(file = here::here("R", "loadLibrary.R"))
source(file = here::here("R", "gc_fake_data.R"))
wakefield::table_heat(x = fakedata, palette = "Set1", flip = TRUE, print = TRUE)
library(readxl)
mydata <- readxl::read_excel(here::here("data", "mydata.xlsx"))
# View(mydata) # Use to view data after importing
# https://cran.r-project.org/web/packages/rio/vignettes/rio.html
# rio::install_formats()

x <- rio::import("mtcars.csv")
y <- rio::import("mtcars.rds")
z <- rio::import("mtcars.dta")

rio::import("mtcars_noext", format = "csv")

rio::export(mtcars, "mtcars.csv")
rio::export(mtcars, "mtcars.rds")
rio::export(mtcars, "mtcars.dta")

rio::export(list(mtcars = mtcars, iris = iris), "multi.xlsx")

# Dataframe report
mydata %>% 
  dplyr::select(-contains("Date")) %>%
  report::report(.)
mydata %>% explore::describe_tbl()
dput(names(mydata))
keycolumns <-  
    mydata %>%  
    sapply(., FUN = dataMaid::isKey) %>%  
    tibble::as_tibble() %>%  
    dplyr::select(  
        which(.[1, ] == TRUE)  
    ) %>%   
    names()  
keycolumns  
mydata %>% 
  dplyr::select(-keycolumns) %>% 
inspectdf::inspect_types()
mydata %>% 
    dplyr::select(-keycolumns,
           -contains("Date")) %>% 
  describer::describe() %>% 
  knitr::kable(format = "markdown")
mydata %>% 
    dplyr::select(-keycolumns) %>% 
  inspectdf::inspect_types() %>% 
  inspectdf::show_plot()
# https://github.com/ropensci/visdat
# http://visdat.njtierney.com/articles/using_visdat.html
# https://cran.r-project.org/web/packages/visdat/index.html
# http://visdat.njtierney.com/

# visdat::vis_guess(mydata)

visdat::vis_dat(mydata)
mydata %>% explore::explore_tbl()
mydata %>% 
    dplyr::select(-keycolumns) %>% 
    inspectdf::inspect_types() %>% 
    dplyr::filter(type == "character") %>% 
    dplyr::select(col_name) %>% 
    dplyr::pull() %>% 
    unlist() -> characterVariables

characterVariables
mydata %>%
    dplyr::select(-keycolumns,
                  -contains("Date")
                  ) %>%
  describer::describe() %>% 
    janitor::clean_names() %>% 
    dplyr::filter(column_type == "factor") %>% 
    dplyr::select(column_name) %>% 
    dplyr::pull() -> categoricalVariables

categoricalVariables
mydata %>%
    dplyr::select(-keycolumns,
                  -contains("Date")) %>%
  describer::describe() %>% 
    janitor::clean_names() %>% 
    dplyr::filter(column_type == "numeric" | column_type == "double") %>% 
    dplyr::select(column_name) %>% 
    dplyr::pull() -> continiousVariables

continiousVariables
mydata %>% 
    dplyr::select(-keycolumns) %>% 
inspectdf::inspect_types() %>% 
  dplyr::filter(type == "numeric") %>% 
  dplyr::select(col_name) %>% 
  dplyr::pull() %>% 
  unlist() -> numericVariables

numericVariables
mydata %>% 
    dplyr::select(-keycolumns) %>% 
inspectdf::inspect_types() %>% 
  dplyr::filter(type == "integer") %>% 
  dplyr::select(col_name) %>% 
  dplyr::pull() %>% 
  unlist() -> integerVariables

integerVariables
mydata %>% 
    dplyr::select(-keycolumns) %>% 
inspectdf::inspect_types() %>% 
  dplyr::filter(type == "list") %>% 
  dplyr::select(col_name) %>% 
  dplyr::pull() %>% 
  unlist() -> listVariables
listVariables
is_date <- function(x) inherits(x, c("POSIXct", "POSIXt"))

dateVariables <- 
names(which(sapply(mydata, FUN = is_date) == TRUE))
dateVariables
View(mydata)
reactable::reactable(data = mydata, sortable = TRUE, resizable = TRUE, filterable = TRUE, searchable = TRUE, pagination = TRUE, paginationType = "numbers", showPageSizeOptions = TRUE, highlight = TRUE, striped = TRUE, outlined = TRUE, compact = TRUE, wrap = FALSE, showSortIcon = TRUE, showSortable = TRUE)
summarytools::view(summarytools::dfSummary(mydata %>% dplyr::select(-keycolumns)))
if(!dir.exists(here::here("out"))) {dir.create(here::here("out"))}

summarytools::view(
  x = summarytools::dfSummary(
    mydata %>% 
      dplyr::select(-keycolumns)
    ),
  file = here::here("out", "mydata_summary.html")
)
if(!dir.exists(here::here("out"))) {dir.create(here::here("out"))}

dataMaid::makeDataReport(data = mydata, 
                         file = here::here("out", "dataMaid_mydata.Rmd"),
                         replace = TRUE,
                         openResult = FALSE, 
                         render = FALSE,
                         quiet = TRUE
                         )
if(!dir.exists(here::here("out"))) {dir.create(here::here("out"))}

mydata %>% 
  dplyr::select(
    -dateVariables
  ) %>% 
  explore::report(
    output_file = "mydata_report.html",
    output_dir = here::here("out") 
    )
dplyr::glimpse(mydata %>% dplyr::select(-keycolumns, -dateVariables))
mydata %>% explore::describe()
explore::explore(mydata)
mydata %>%
  explore::explore_all()
visdat::vis_expect(data = mydata,
                   expectation = ~.x == -1,
                   show_perc = TRUE)

visdat::vis_expect(mydata, ~.x >= 25)
visdat::vis_miss(airquality,
                 cluster = TRUE)
visdat::vis_miss(airquality,
         sort_miss = TRUE)
xray::anomalies(mydata)
xray::distributions(mydata)
DataExplorer::plot_str(mydata)
DataExplorer::plot_str(mydata, type = "r")
DataExplorer::introduce(mydata)
DataExplorer::plot_intro(mydata)
DataExplorer::plot_missing(mydata)
mydata2 <- DataExplorer::drop_columns(mydata, "TStage")
DataExplorer::plot_bar(mydata)
DataExplorer::plot_bar(mydata, with = "Death")
DataExplorer::plot_histogram(mydata)
if(!dir.exists(here::here("out"))) {dir.create(here::here("out"))}

# https://cran.r-project.org/web/packages/dataMaid/vignettes/extending_dataMaid.html
library("dataMaid")
dataMaid::makeDataReport(mydata,
  #add extra precheck function
  preChecks = c("isKey", "isSingular", "isSupported", "isID"),

  #Add the extra summaries - countZeros() for character, factor,
  #integer, labelled and numeric variables and meanSummary() for integer,
  #numeric and logical variables:
  summaries = setSummaries(
    character = defaultCharacterSummaries(add = "countZeros"),
    factor = defaultFactorSummaries(add = "countZeros"),
    labelled = defaultLabelledSummaries(add = "countZeros"),
    numeric = defaultNumericSummaries(add = c("countZeros", "meanSummary")),
    integer = defaultIntegerSummaries(add = c("countZeros", "meanSummary")),
    logical = defaultLogicalSummaries(add =  c("meanSummary"))
  ),

  #choose mosaicVisual() for categorical variables,
  #prettierHist() for all others:
  visuals = setVisuals(
    factor = "mosaicVisual",
    numeric = "prettierHist",
    integer = "prettierHist",
    Date = "prettierHist"
  ),

  #Add the new checkFunction, identifyColons, for character, factor and
  #labelled variables:
  checks = setChecks(
    character = defaultCharacterChecks(add = "identifyColons"),
    factor = defaultFactorChecks(add = "identifyColons"),
    labelled = defaultLabelledChecks(add = "identifyColons")
  ),

  #overwrite old versions of the report, render to html and don't
  #open the html file automatically:
  replace = TRUE,
  output = "html",
  open = FALSE,
  file = here::here("out/dataMaid_mydata.Rmd")
)
# https://cran.r-project.org/web/packages/summarytools/vignettes/Recommendations-rmarkdown.html
# https://github.com/dcomtois/summarytools
library(knitr)
opts_chunk$set(comment=NA,
               prompt=FALSE,
               cache=FALSE,
               echo=TRUE,
               results='asis' # add to individual summarytools chunks
               )
library(summarytools)
st_css()
st_options(bootstrap.css     = FALSE,       # Already part of the theme so no need for it
           plain.ascii       = FALSE,       # One of the essential settings
           style             = "rmarkdown", # Idem.
           dfSummary.silent  = TRUE,        # Suppresses messages about temporary files
           footnote          = NA,          # Keeping the results minimalistic
           subtitle.emphasis = FALSE)       # For the vignette theme, this gives
                                            # much better results. Your mileage may vary.
summarytools::freq(iris$Species, plain.ascii = FALSE, style = "rmarkdown")
summarytools::freq(iris$Species, report.nas = FALSE, headings = FALSE, cumul = TRUE, totals = TRUE)
summarytools::freq(tobacco$gender, style = 'rmarkdown')
summarytools::freq(tobacco[ ,c("gender", "age.gr", "smoker")])
print(freq(tobacco$gender), method = 'render')
view(dfSummary(iris))
dfSummary(tobacco, style = 'grid', graph.magnif = 0.75, tmp.img.dir = "/tmp")
dfSummary(tobacco, plain.ascii = FALSE, style = "grid",
          graph.magnif = 0.75, valid.col = FALSE, tmp.img.dir = "/tmp")
print(dfSummary(tobacco, graph.magnif = 0.75), method = 'render')
# https://github.com/rolkra/explore
# https://cran.r-project.org/web/packages/explore/vignettes/explore.html
# https://cran.r-project.org/web/packages/explore/vignettes/explore_mtcars.html


# library(dplyr)
# library(explore)

explore::explore(mydata)

# iris %>% report(output_file = "report.html", output_dir = here::here())


# iris$is_versicolor <- ifelse(iris$Species == "versicolor", 1, 0)
# iris %>%
#   report(output_file = "report.html",
#          output_dir = here::here(),
#          target = is_versicolor
# # , split = FALSE
# )
iris %>% explore::explore_tbl()

iris %>% explore::describe_tbl()

iris %>% explore::explore(Species)
iris %>% explore::explore(Sepal.Length)
iris %>% explore::explore(Sepal.Length, target = is_versicolor)
iris %>% explore::explore(Sepal.Length, target = is_versicolor, split = FALSE)
iris %>% explore::explore(Sepal.Length, target = Species)
iris %>% explore::explore(Sepal.Length, target = Petal.Length)


%>% %>%
  explore::explore_all()

iris %>%
  dplyr::select(Sepal.Length, Sepal.Width) %>%
  explore::explore_all()

iris %>%
  dplyr::select(Sepal.Length, Sepal.Width, is_versicolor) %>%
  explore::explore_all(target = is_versicolor)

iris %>%
  dplyr::select(Sepal.Length, Sepal.Width, is_versicolor) %>%
  explore::explore_all(target = is_versicolor, split = FALSE)

iris %>%
  dplyr::select(Sepal.Length, Sepal.Width, Species) %>%
  explore::explore_all(target = Species)

iris %>%
  dplyr::select(Sepal.Length, Sepal.Width, Petal.Length) %>%
  explore::explore_all(target = Petal.Length)
iris %>%
  explore::explore_all()

knitr::opts_current(fig.height=explore::total_fig_height(iris, target = Species))

explore::total_fig_height(iris, target = Species)

iris %>% explore::explore_all(target = Species)
iris %>% explore::explore(Sepal.Length, min_val = 4.5, max_val = 7)
iris %>% explore::explore(Sepal.Length, auto_scale = FALSE)
mtcars %>% explore::describe()
# https://cran.r-project.org/web/packages/dlookr/vignettes/EDA.html

dlookr::describe(mydata
                 # ,
                 # cols = c(statistic)
                 )

# dlookr::describe(carseats, Sales, CompPrice, Income)
# dlookr::describe(carseats, Sales:Income)
# dlookr::describe(carseats, -(Sales:Income))
mydata %>%
  dlookr::describe() %>%
  dplyr::select(variable, skewness, mean, p25, p50, p75) %>%
  dplyr::filter(!is.na(skewness)) %>%
  arrange(desc(abs(skewness)))
# https://cran.r-project.org/web/packages/dlookr/vignettes/EDA.html

carseats %>%
  dlookr::eda_report(target = Sales,
                     output_format = "pdf",
             output_file = "EDA.pdf"
             )
carseats %>%
  dlookr::eda_report(target = Sales,
             output_format = "html",
             output_file = "EDA.html"
             )
# install.packages("ISLR")
library("ISLR")
# install.packages("SmartEDA")
library("SmartEDA")
## Load sample dataset from ISLR pacakge
Carseats <- ISLR::Carseats

## overview of the data;
SmartEDA::ExpData(data=Carseats,type=1)
## structure of the data
SmartEDA::ExpData(data=Carseats,type=2)


mydata %>%
  dplyr::select(-dplyr::contains("Date")) %>%
  report::report()
# cat(names(mydata), sep = " + \n")
library(arsenal)
tab1 <- arsenal::tableby(
  ~ Sex +
    Age +
    Race +
    PreinvasiveComponent +
    LVI +
    PNI +
    Death +
    Group +
    Grade +
    TStage +
    # `Anti-X-intensity` +
    # `Anti-Y-intensity` +
    LymphNodeMetastasis +
    Valid +
    Smoker +
    Grade_Level
  ,
  data = mydata 
)
summary(tab1)
library(tableone)
mydata %>% 
  dplyr::select(-keycolumns,
         -dateVariables
         ) %>% 
tableone::CreateTableOne(data = .)
# CreateTableOne(vars = myVars, data = mydata, factorVars = characterVariables)
# tab <- CreateTableOne(vars = myVars, data = pbc, factorVars = catVars)
# print(tab, showAllLevels = TRUE)
# ?print.TableOne
# summary(tab)
# print(tab, nonnormal = biomarkers)
# print(tab, nonnormal = biomarkers, exact = "stage", quote = TRUE, noSpaces = TRUE)
# tab3Mat <- print(tab3, nonnormal = biomarkers, exact = "stage", quote = FALSE, noSpaces = TRUE, printToggle = FALSE)
# write.csv(tab3Mat, file = "myTable.csv")
mydata %>% 
  dplyr::select(
    continiousVariables,
    numericVariables,
    integerVariables
  ) %>% 
summarytools::descr(., style = 'rmarkdown')
print(summarytools::descr(mydata), method = 'render', table.classes = 'st-small')
mydata %>% 
  summarytools::descr(.,
                      stats = "common",
                      transpose = TRUE,
                      headings = FALSE
                      )
mydata %>% 
  summarytools::descr(stats = "common") %>%
  summarytools::tb()
mydata$Sex %>% 
  summarytools::freq(cumul = FALSE, report.nas = FALSE) %>%
  summarytools::tb()
mydata %>%
  explore::describe() %>%
  dplyr::filter(unique < 5)
mydata %>%
  explore::describe() %>%
  dplyr::filter(na > 0)
mydata %>% explore::describe()
source(here::here("R", "gc_desc_cat.R"))
tab <- 
  mydata %>% 
  dplyr::select(
    -keycolumns
    ) %>% 
  tableone::CreateTableOne(data = .)
?print.CatTable
tab$CatTable
race_stats <- summarytools::freq(mydata$Race) 
print(race_stats,
      report.nas = FALSE,
      totals = FALSE,
      display.type = FALSE,
      Variable.label = "Race Group"
      )
mydata %>% explore::describe(PreinvasiveComponent)
## Frequency or custom tables for categorical variables
SmartEDA::ExpCTable(
  mydata,
  Target = NULL,
  margin = 1,
  clim = 10,
  nlim = 5,
  round = 2,
  bin = NULL,
  per = T
)
inspectdf::inspect_cat(mydata)

inspectdf::inspect_cat(mydata)$levels$Group
library(summarytools)

grouped_freqs <- stby(data = mydata$Smoker, 
                      INDICES = mydata$Sex, 
                      FUN = freq, cumul = FALSE, report.nas = FALSE)

grouped_freqs %>% tb(order = 2)
summarytools::stby(
  list(x = mydata$LVI, y = mydata$LymphNodeMetastasis), 
  mydata$PNI,
  summarytools::ctable
  )
with(mydata, 
     summarytools::stby(
       list(x = LVI, y = LymphNodeMetastasis), PNI,
       summarytools::ctable
       )
     )
SmartEDA::ExpCTable(
  mydata,
  Target = "Sex",
  margin = 1,
  clim = 10,
  nlim = NULL,
  round = 2,
  bin = 4,
  per = F
)
mydata %>% 
  dplyr::select(characterVariables) %>% 
  dplyr::select(PreinvasiveComponent,
         PNI,
         LVI
         ) %>% 
reactable::reactable(data = ., groupBy = c("PreinvasiveComponent", "PNI"), columns = list(
  LVI = reactable::colDef(aggregate = "count")
))
questionr:::icut()
source(here::here("R", "gc_desc_cont.R"))
tab <- tableone::CreateTableOne(data = mydata)
# ?print.ContTable
tab$ContTable
print(tab$ContTable, nonnormal = c("Anti-X-intensity"))
mydata %>% explore::describe(Age)
mydata %>% 
  dplyr::select(continiousVariables) %>% 
SmartEDA::ExpNumStat(
  data = .,
  by = "A",
  gp = NULL,
  Qnt = seq(0, 1, 0.1),
  MesofShape = 2,
  Outlier = TRUE,
  round = 2
)
inspectdf::inspect_num(mydata, breaks = 10)
inspectdf::inspect_num(mydata)$hist$Age
inspectdf::inspect_num(mydata, breaks = 10) %>%
  inspectdf::show_plot()
grouped_descr <- summarytools::stby(data = mydata, 
                      INDICES = mydata$Sex, 
                      FUN = summarytools::descr, stats = "common")
# grouped_descr %>% summarytools::tb(order = 2)
grouped_descr %>% summarytools::tb()
mydata %>%
  group_by(US) %>% 
  dlookr::describe(Sales, Income) 

carseats %>%
  group_by(US, Urban) %>% 
  dlookr::describe(Sales, Income) 

categ <- dlookr::target_by(carseats, US)
cat_num <- dlookr::relate(categ, Sales)
cat_num
summary(cat_num)
plot(cat_num)
   summarytools::stby(data = mydata, 
                               INDICES = mydata$PreinvasiveComponent, 
                               FUN = summarytools::descr,
                      stats = c("mean", "sd", "min", "med", "max"), 
                               transpose = TRUE)
with(mydata, 
     summarytools::stby(Age, PreinvasiveComponent, summarytools::descr), 
                   stats = c("mean", "sd", "min", "med", "max"),
                   transpose = TRUE
                   )
mydata %>% 
  group_by(PreinvasiveComponent) %>% 
  summarytools::descr(stats = "fivenum")
## Summary statistics by – category
SmartEDA::ExpNumStat(
  mydata,
  by = "GA",
  gp = "PreinvasiveComponent",
  Qnt = seq(0, 1, 0.1),
  MesofShape = 2,
  Outlier = TRUE,
  round = 2
)
mydata %>% 
  janitor::tabyl(Sex) %>%
  janitor::adorn_pct_formatting(rounding = 'half up', digits = 1) %>%
  knitr::kable()
mydata %>% 
  janitor::tabyl(Race) %>%
  janitor::adorn_pct_formatting(rounding = 'half up', digits = 1) %>%
  knitr::kable()
mydata %>% 
  janitor::tabyl(PreinvasiveComponent) %>%
  janitor::adorn_pct_formatting(rounding = 'half up', digits = 1) %>%
  knitr::kable()
mydata %>% 
  janitor::tabyl(LVI) %>%
  janitor::adorn_pct_formatting(rounding = 'half up', digits = 1) %>%
  knitr::kable()
mydata %>% 
  janitor::tabyl(PNI) %>%
  janitor::adorn_pct_formatting(rounding = 'half up', digits = 1) %>%
  knitr::kable()
mydata %>% 
  janitor::tabyl(Group) %>%
  janitor::adorn_pct_formatting(rounding = 'half up', digits = 1) %>%
  knitr::kable()
mydata %>% 
  janitor::tabyl(Grade) %>%
  janitor::adorn_pct_formatting(rounding = 'half up', digits = 1) %>%
  knitr::kable()
mydata %>% 
  janitor::tabyl(TStage) %>%
  janitor::adorn_pct_formatting(rounding = 'half up', digits = 1) %>%
  knitr::kable()
mydata %>% 
  janitor::tabyl(LymphNodeMetastasis) %>%
  janitor::adorn_pct_formatting(rounding = 'half up', digits = 1) %>%
  knitr::kable()
mydata %>% 
  janitor::tabyl(Grade_Level) %>%
  janitor::adorn_pct_formatting(rounding = 'half up', digits = 1) %>%
  knitr::kable()
mydata %>% 
  janitor::tabyl(DeathTime) %>%
  janitor::adorn_pct_formatting(rounding = 'half up', digits = 1) %>%
  knitr::kable()
mydata %>% 
jmv::descriptives(
    data = .,
    vars = 'Age',
    hist = TRUE,
    dens = TRUE,
    box = TRUE,
    violin = TRUE,
    dot = TRUE,
    mode = TRUE,
    sd = TRUE,
    variance = TRUE,
    skew = TRUE,
    kurt = TRUE,
    quart = TRUE)
mydata %>% 
jmv::descriptives(
    data = .,
    vars = 'AntiX_intensity',
    hist = TRUE,
    dens = TRUE,
    box = TRUE,
    violin = TRUE,
    dot = TRUE,
    mode = TRUE,
    sd = TRUE,
    variance = TRUE,
    skew = TRUE,
    kurt = TRUE,
    quart = TRUE)
mydata %>% 
jmv::descriptives(
    data = .,
    vars = 'AntiY_intensity',
    hist = TRUE,
    dens = TRUE,
    box = TRUE,
    violin = TRUE,
    dot = TRUE,
    mode = TRUE,
    sd = TRUE,
    variance = TRUE,
    skew = TRUE,
    kurt = TRUE,
    quart = TRUE)
library(finalfit)
dependent <- c("dependent1",
               "dependent2"
              )

explanatory <- c("explanatory1",
                 "explanatory2"
                 )

dependent <- "PreinvasiveComponent"

explanatory <- c("Sex", "Age", "Grade", "TStage")

source(here::here("R", "gc_table_cross.R"))
CreateTableOne(vars = myVars, strata = "columnname", data = pbc, factorVars = catVars)
print(tab, nonnormal = biomarkers, exact = "exactVariable", smd = TRUE)

write2html(
  knitr::kable(head(mockstudy)), paste0(tmpdir, "/test.kable.keep.rmd.html"),
  quiet = TRUE, # passed to rmarkdown::render
  keep.rmd = TRUE
)
ctable(tobacco$gender, tobacco$smoker, style = 'rmarkdown')
print(ctable(tobacco$gender, tobacco$smoker), method = 'render')
print(ctable(tobacco$smoker, tobacco$diseased, prop = "r"), method = "render")
with(tobacco, 
     print(ctable(smoker, diseased, prop = 'n', totals = FALSE, chisq = TRUE),
           headings = FALSE, method = "render"))
# devtools::install_github("ewenharrison/summarizer")
# library(summarizer)
# data(colon_s)
explanatory = c("age", "age.factor", "sex.factor", "obstruct.factor")
dependent = "perfor.factor"
colon_s %>%
  summary.factorlist(dependent, explanatory, p=TRUE) %>% 
    knitr::kable(row.names=FALSE, align=c("l", "l", "r", "r", "r"))

explanatory = c("age.factor", "sex.factor", "obstruct.factor", "perfor.factor")
dependent = 'mort_5yr'
colon_s %>%
  summary.factorlist(dependent, explanatory) %>% 
    knitr::kable(row.names=FALSE, align=c("l", "l", "r", "r", "r"))


mydata %>%
    summary_factorlist(dependent = 'PreinvasiveComponent', 
                       explanatory = explanatory,
                       # column = TRUE,
                       total_col = TRUE,
                       p = TRUE,
                       add_dependent_label = TRUE,
                       na_include=FALSE
                       # catTest = catTestfisher
                       ) -> table

knitr::kable(table, row.names = FALSE, align = c('l', 'l', 'r', 'r', 'r'))

table1 <- arsenal::tableby(PreinvasiveComponent ~ explanatory, mydata)

summary(table1)

knitr::kable(table1,
                         row.names = FALSE,
                         align = c('l', 'l', 'r', 'r', 'r', 'r'),
                         format = 'html') %>%
                kableExtra::kable_styling(kable_input = .,
                                          bootstrap_options = 'striped',
                                          full_width = F,
                                          position = 'left')
tangram::tangram(PreinvasiveComponent ~ explanatory, mydata)

tangram::html5(tangram::tangram(PreinvasiveComponent ~ explanatory, mydata),
                    fragment = TRUE,
                    inline = 'nejm.css',
                    caption = 'Cross TablePreinvasiveComponentNEJM Style',
                    id = 'tbl3')
tangram::html5(tangram::tangram(PreinvasiveComponent ~ explanatory, mydata),
                    fragment = TRUE,
                    inline = 'lancet.css',
                    caption = 'Cross TablePreinvasiveComponentLancet Style',
                    id = 'tbl3')
mydata$int <- lubridate::interval(
  lubridate::ymd(mydata$SurgeryDate),
  lubridate::ymd(mydata$LastFollowUpDate)
  )
mydata$OverallTime <- lubridate::time_length(mydata$int, "month")
mydata$OverallTime <- round(mydata$OverallTime, digits = 1)
mydata$OverallTime <- mydata$genel_sagkalim
## Recoding mydata$Death into mydata$Outcome
mydata$Outcome <- forcats::fct_recode(as.character(mydata$Death),
               "1" = "TRUE",
               "0" = "FALSE")
mydata$Outcome <- as.numeric(as.character(mydata$Outcome))
table(mydata$Death, mydata$Outcome)
library(survival)
# data(lung)
# km <- with(lung, Surv(time, status))
km <- with(mydata, Surv(OverallTime, Outcome))
head(km,80)
plot(km)
# Drawing Survival Curves Using ggplot2
# https://rpkgs.datanovia.com/survminer/reference/ggsurvplot.html
dependentKM <- "Surv(OverallTime, Outcome)"
explanatoryKM <- "LVI"

mydata %>%
  finalfit::surv_plot(.data = .,
                      dependent = dependentKM,
                      explanatory = explanatoryKM,
                      xlab='Time (months)',
                      pval=TRUE,
                      legend = 'none',
                      break.time.by = 12,
                      xlim = c(0,60)
                      # legend.labs = c('a','b')
                      )
# Drawing Survival Curves Using ggplot2
# https://rpkgs.datanovia.com/survminer/reference/ggsurvplot.html

mydata %>%
  finalfit::surv_plot(.data = .,
                      dependent = "Surv(OverallTime, Outcome)",
                      explanatory = "LVI",
                      xlab='Time (months)',
                      pval=TRUE,
                      legend = 'none',
                      break.time.by = 12,
                      xlim = c(0,60)
                      # legend.labs = c('a','b')
                      )
library(finalfit)
library(survival)
explanatoryUni <- "LVI"
dependentUni <- "Surv(OverallTime, Outcome)"

mydata %>%
finalfit::finalfit(dependentUni, explanatoryUni) -> tUni

knitr::kable(tUni, row.names=FALSE, align=c('l', 'l', 'r', 'r', 'r', 'r'))
tUni_df <- tibble::as_tibble(tUni, .name_repair = "minimal") %>% 
  janitor::clean_names() 

tUni_df_descr <- paste0("When ",
                        tUni_df$dependent_surv_overall_time_outcome[1],
                        " is ",
                        tUni_df$x[2],
                        ", there is ",
                        tUni_df$hr_univariable[2],
                        " times risk than ",
                        "when ",
                        tUni_df$dependent_surv_overall_time_outcome[1],
                        " is ",
                        tUni_df$x[1],
                        "."
                        )
km_fit <- survfit(Surv(OverallTime, Outcome) ~ LVI, data = mydata)
km_fit
plot(km_fit)
# summary(km_fit)
km_fit_median_df <- summary(km_fit)
km_fit_median_df <- as.data.frame(km_fit_median_df$table) %>% 
  janitor::clean_names() %>% 
  tibble::rownames_to_column()
km_fit_median_df <- summary(km_fit)
            km_fit_median_df <- as.data.frame(km_fit_median_df$table) %>%
                tibble::rownames_to_column()
            
            names(km_fit_median_df) <- paste0("m", 1:dim(km_fit_median_df)[2])

            km_fit_median_definition <- 
                
            km_fit_median_df %>%
                dplyr::mutate(
                    description =
                        glue::glue(
                            "When {m1}, median survival is {m8} [{m9} - {m10}, 95% CI] months."
                        )
                ) %>%
                dplyr::select(description) %>%
                dplyr::pull() 
sTable <- summary(km_fit)$table
            st <- data.frame()
            
            for (i in seq_len(nrow(km_fit))) {
                if (nrow(km_fit) == 1)
                    g <- sTable
                else
                    g <- sTable[i,]
                nevents <- sum(g['events'])
                n <- g['n.max']
                ncensor <- n - nevents
                median <- g['median']
                mean <- g['*rmean']
                prop <- nevents / n
                
                print(rowNo=i, list(
                    censored=ncensor,
                    events=nevents,
                    n=n,
                    prop=nevents/n,
                    median=median,
                    mean=mean))
            }
            
            st$setStatus('complete')
            
            
            results1 <- st
km_fit
broom::tidy(km_fit)
km_fit_median_df %>% 
  dplyr::mutate(
    description = 
      glue::glue(
      "When {rowname}, median survival is {median} [{x0_95lcl} - {x0_95ucl}, 95% CI] months."
    )
  ) %>% 
  dplyr::select(description) %>% 
  dplyr::pull() -> km_fit_median_definition
summary(km_fit, times = c(12,36,60))
km_fit_summary <- summary(km_fit, times = c(12,36,60))

km_fit_df <- as.data.frame(km_fit_summary[c("strata", "time", "n.risk", "n.event", "surv", "std.err", "lower", "upper")])
km_fit_df %>% 
  dplyr::mutate(
    description = 
      glue::glue(
      "When {strata}, {time} month survival is {scales::percent(surv)} [{scales::percent(lower)}-{scales::percent(upper)}, 95% CI]."
    )
  ) %>% 
  dplyr::select(description) %>% 
  dplyr::pull() -> km_fit_definition
library(survival)
surv_fit <- survival::survfit(Surv(time, status) ~ ph.ecog, data=lung)
insight::is_model_supported(surv_fit)
insight::find_formula(surv_fit)
report::report_participants(mydata)
dependentKM <-  "Surv(OverallTime, Outcome2)"

explanatoryKM <- c("explanatory1",
               "explanatory2"
               )
source(here::here("R", "gc_survival.R"))
mydependent <-  "Surv(time, status)"
explanatory <- "Organ"

mysurvival <- function(mydata, mydependent, explanatory) {
    {{mydata}} %>%
        finalfit::surv_plot(dependent = {{mydependent}},
                            explanatory = {{explanatory}},
                            xlab='Time (months)',
                            pval=TRUE,
                            legend = 'none',
                            break.time.by = 12,
                            xlim = c(0,60)
        )
}


# library(tidyverse)
mysurvival(mydata = whippleables, mydependent = mydependent, explanatory = explanatory)

explanatory <- c("Organ", "LVI")

deneme <- purrr::map(explanatory, mysurvival, mydata = whippleables, mydependent = mydependent)

dependentKM <- "Surv(OverallTime, Outcome)"
explanatoryKM <- "TStage"

mydata %>%
  finalfit::surv_plot(.data = .,
                      dependent = dependentKM,
                      explanatory = explanatoryKM,
                      xlab='Time (months)',
                      pval=TRUE,
                      legend = 'none',
                      break.time.by = 12,
                      xlim = c(0,60)
                      # legend.labs = c('a','b')
                      )
survminer::pairwise_survdiff(
  formula = Surv(OverallTime, Outcome) ~ TStage, 
                             data = mydata,
                             p.adjust.method = "BH"
  )
km_fit
print(km_fit, 
      scale=1,
      digits = max(options()$digits - 4,3),
      print.rmean=getOption("survfit.print.rmean"),
      rmean = getOption('survfit.rmean'),
      print.median=getOption("survfit.print.median"),
      median = getOption('survfit.median')
      
      )
library(finalfit)
library(survival)
explanatoryMultivariate <- explanatoryKM
dependentMultivariate <- dependentKM

mydata %>%
  finalfit(dependentMultivariate, explanatoryMultivariate) -> tMultivariate

knitr::kable(tMultivariate, row.names=FALSE, align=c("l", "l", "r", "r", "r", "r"))
# https://tidymodels.github.io/parsnip/reference/surv_reg.html
library(parsnip)
surv_reg()
#> Parametric Survival Regression Model Specification (regression)
#> # Parameters can be represented by a placeholder:
surv_reg(dist = varying())

#> Parametric Survival Regression Model Specification (regression)
#> 
#> Main Arguments:
#>   dist = varying()
#> 
model <- surv_reg(dist = "weibull")
model
#> Parametric Survival Regression Model Specification (regression)
#> 
#> Main Arguments:
#>   dist = weibull
#> update(model, dist = "lnorm")#> Parametric Survival Regression Model Specification (regression)
#> 
#> Main Arguments:
#>   dist = lnorm
#> 


# From randomForest
rf_1 <- randomForest(x, y, mtry = 12, ntree = 2000, importance = TRUE)

# From ranger
rf_2 <- ranger(
  y ~ ., 
  data = dat, 
  mtry = 12, 
  num.trees = 2000, 
  importance = 'impurity'
)

# From sparklyr
rf_3 <- ml_random_forest(
  dat, 
  intercept = FALSE, 
  response = "y", 
  features = names(dat)[names(dat) != "y"], 
  col.sample.rate = 12,
  num.trees = 2000
)




rand_forest(mtry = 12, trees = 2000) %>%
  set_engine("ranger", importance = 'impurity') %>%
  fit(y ~ ., data = dat)


rand_forest(mtry = 12, trees = 2000) %>%
  set_engine("spark") %>%
  fit(y ~ ., data = dat)



library(survival)
library(survminer)
library(finalfit)

mydata %>%
  finalfit::surv_plot('Surv(OverallTime, Outcome)', 'LVI', 
  xlab='Time (months)', pval=TRUE, legend = 'none',
    break.time.by = 12, xlim = c(0,60)

# legend.labs = c('a','b')

)
explanatoryUni <- 'LVI'
dependentUni <- 'Surv(OverallTime, Outcome)'
mydata %>%
finalfit(dependentUni, explanatoryUni) -> tUni

knitr::kable(tUni, row.names=FALSE, align=c('l', 'l', 'r', 'r', 'r', 'r'))

tUni_df <- tibble::as_tibble(tUni, .name_repair = 'minimal') %>%
janitor::clean_names(dat = ., case = 'snake')


n_level <- dim(tUni_df)[1]

tUni_df_descr <- function(n) {
    paste0(
        'When ',
        tUni_df$dependent_surv_overall_time_outcome[1],
        ' is ',
        tUni_df$x[n + 1],
        ', there is ',
        tUni_df$hr_univariable[n + 1],
        ' times risk than ',
        'when ',
        tUni_df$dependent_surv_overall_time_outcome[1],
        ' is ',
        tUni_df$x[1],
        '.'
    )

}



results5 <- purrr::map(.x = c(2:n_level-1), .f = tUni_df_descr)

print(unlist(results5))

km_fit <- survfit(Surv(OverallTime, Outcome) ~ LVI, data = mydata)
km_fit

km_fit_median_df <- summary(km_fit)
km_fit_median_df <- as.data.frame(km_fit_median_df$table) %>%
    janitor::clean_names(dat = ., case = 'snake') %>%
    tibble::rownames_to_column(.data = ., var = 'LVI')



km_fit_median_df %>%
    dplyr::mutate(
        description =
        glue::glue(
        'When, LVI, {LVI}, median survival is {median} [{x0_95lcl} - {x0_95ucl}, 95% CI] months.'
)
    ) %>%
        dplyr::mutate(
description = gsub(pattern = 'thefactor=', replacement = ' is ', x = description)
        ) %>%
    dplyr::select(description) %>%
    dplyr::pull() -> km_fit_median_definition

km_fit_median_definition




summary(km_fit, times = c(12,36,60))

km_fit_summary <- summary(km_fit, times = c(12,36,60))

km_fit_df <- as.data.frame(km_fit_summary[c('strata', 'time', 'n.risk', 'n.event', 'surv', 'std.err', 'lower', 'upper')])

km_fit_df

     


km_fit_df %>%
    dplyr::mutate(
        description =
glue::glue(
    'When {strata}, {time} month survival is {scales::percent(surv)} [{scales::percent(lower)}-{scales::percent(upper)}, 95% CI].'
)
    ) %>%
    dplyr::select(description) %>%
    dplyr::pull() -> km_fit_definition

km_fit_definition


summary(km_fit)$table

km_fit_median_df <- summary(km_fit)
results1html <- as.data.frame(km_fit_median_df$table) %>%
    janitor::clean_names(dat = ., case = 'snake') %>%
    tibble::rownames_to_column(.data = ., var = 'LVI')

results1html[,1] <- gsub(pattern = 'thefactor=',
 replacement = '',
 x = results1html[,1])

knitr::kable(results1html,
 row.names = FALSE,
 align = c('l', rep('r', 9)),
 format = 'html',
 digits = 1)

     

formula_p1 <- jmvcore::constructFormula(terms = self$options$overalltime)

formula_p3 <- jmvcore::constructFormula(terms = LVI)

formula_p2 <- jmvcore::constructFormula(terms = self$options$outcome)


formula_p <- paste('Surv(', formula_p1, ',',  formula_p2, ') ~ ', formula_p3)

formula_p <- as.formula(formula_p)


results8 <-
    survminer::pairwise_survdiff(
    formula = formula_p,
    data = self$data,
    p.adjust.method = 'BH'
)
projectName <- list.files(path = here::here(), pattern = "Rproj")
projectName <- gsub(pattern = ".Rproj", replacement = "", x = projectName)

analysisDate <- as.character(Sys.Date())

imageName <- paste0(projectName, analysisDate, ".RData")

save.image(file = here::here("data", imageName))

rdsName <- paste0(projectName, analysisDate, ".rds")

readr::write_rds(x = mydata, path = here::here("data", rdsName))

saveRDS(object = mydata, file = here::here("data", rdsName))

excelName <- paste0(projectName, analysisDate, ".xlsx")

rio::export(
  x = mydata,
  file = here::here("data", excelName),
  format = "xlsx"
)

# writexl::write_xlsx(mydata, here::here("data", excelName))

print(glue::glue(
    "saved data after analysis to ",
    rownames(file.info(here::here("data", excelName))),
    " : ",
    as.character(
        file.info(here::here("data", excelName))$ctime
    )
    )
)
mydata %>% 
  downloadthis::download_this(
    output_name = excelName,
    output_extension = ".csv",
    button_label = "Download data as csv",
    button_type = "default"
  )

mydata %>% 
  downloadthis::download_this(
    output_name = excelName,
    output_extension = ".xlsx",
    button_label = "Download data as xlsx",
    button_type = "primary"
  )
# pacman::p_load(here, lubridate, glue)
# here::here("data", glue("{today()}_trends.csv"))
# mydata %>% select(
#     -c(
#         rapor_yil,
#         rapor_no,
#         protokol_no,
#         istek_tarihi,
#         nux_yada_met_varsa_tarihi,
#         son_hastane_vizit_tarihi,
#         Outcome
#     )
# ) -> finalSummary
# 
# summarytools::view(summarytools::dfSummary(x = finalSummary
#                                            , style = "markdown"))
citation()
report::cite_packages(session = sessionInfo())
report::show_packages(session = sessionInfo()) %>% 
    kableExtra::kable()
# citation("tidyverse")
citation("readxl")
citation("janitor")
# citation("report")
citation("finalfit")
# citation("ggstatsplot")
if(!dir.exists(here::here("bib"))) {dir.create(here::here("bib"))}

knitr::write_bib(x = c(.packages(), "knitr", "shiny"),
                 file = here::here("bib", "packages.bib")
)
sessionInfo()
pacman::p_loaded(all = TRUE)
search()
library()
installed.packages()[1:5, c("Package", "Version")]
installed.packages()

---

References

Knijn, N., F. Simmer, and I. D. Nagtegaal. 2015. “Recommendations for Reporting Histopathology Studies: A Proposal.” Virchows Archiv 466 (6): 611–15. https://doi.org/10.1007/s00428-015-1762-3.

Schmidt, Robert L., Deborah J. Chute, Jorie M. Colbert-Getz, Adolfo Firpo-Betancourt, Daniel S. James, Julie K. Karp, Douglas C. Miller, et al. 2017. “Statistical Literacy Among Academic Pathologists: A Survey Study to Gauge Knowledge of Frequently Used Statistical Tests Among Trainees and Faculty.” Archives of Pathology & Laboratory Medicine 141 (2): 279–87. https://doi.org/10.5858/arpa.2016-0200-OA.

---

1. From Table 1: Proposed items for reporting histopathology studies. Recommendations for reporting histopathology studies: a proposal Virchows Arch (2015) 466:611–615 DOI 10.1007/s00428-015-1762-3 https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4460276/↩︎

2. From Table 1: Proposed items for reporting histopathology studies. Recommendations for reporting histopathology studies: a proposal Virchows Arch (2015) 466:611–615 DOI 10.1007/s00428-015-1762-3 https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4460276/↩︎

3. See childRmd/_01header.Rmd file for other general settings↩︎

4. Change echo = FALSE to hide codes after knitting.↩︎

5. See childRmd/_02fakeData.Rmd file for other codes↩︎

6. Synthea The validity of synthetic clinical data: a validation study of a leading synthetic data generator (Synthea) using clinical quality measures. BMC Med Inform Decis Mak 19, 44 (2019) doi:10.1186/s12911-019-0793-0↩︎

7. https://bmcmedinformdecismak.biomedcentral.com/articles/10.1186/s12911-019-0793-0↩︎

8. Synthetic Patient Generation↩︎

9. Basic Setup and Running↩︎

10. intelligent patient data generator (iPDG)↩︎

11. https://medium.com/free-code-camp/how-our-test-data-generator-makes-fake-data-look-real-ace01c5bde4a↩︎

12. https://forums.librehealth.io/t/demo-data-generation/203↩︎

13. https://mihin.org/services/patient-generator/↩︎

14. lung, cancer, breast datası ile birleştir↩︎

15. See childRmd/_03importData.Rmd file for other codes↩︎

16. See childRmd/_04briefSummary.Rmd file for other codes↩︎

17. https://www.hhs.gov/hipaa/index.html↩︎

18. Kişisel verilerin kaydedilmesi ve kişisel verileri hukuka aykırı olarak verme veya ele geçirme Türk Ceza Kanunu’nun 135. ve 136. maddesi kapsamında bizim hukuk sistemimizde suç olarak tanımlanmıştır. Kişisel verilerin kaydedilmesi suçunun cezası 1 ila 3 yıl hapis cezasıdır. Suçun nitelikli hali ise, kamu görevlisi tarafından görevin verdiği yetkinin kötüye kullanılarak veya belirli bir meslek veya sanatın sağladığı kolaylıktan yararlanılarak işlenmesidir ki bu durumda suçun cezası 1.5 ile 4.5 yıl hapis cezası olacaktır.↩︎

19. See childRmd/_06variableTypes.Rmd file for other codes↩︎

20. See childRmd/_07overView.Rmd file for other codes↩︎

21. Statistical Literacy Among Academic Pathologists: A Survey Study to Gauge Knowledge of Frequently Used Statistical Tests Among Trainees and Faculty. Archives of Pathology & Laboratory Medicine: February 2017, Vol. 141, No. 2, pp. 279-287. https://doi.org/10.5858/arpa.2016-0200-OA↩︎

22. From Table 1: Proposed items for reporting histopathology studies. Recommendations for reporting histopathology studies: a proposal Virchows Arch (2015) 466:611–615 DOI 10.1007/s00428-015-1762-3 https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4460276/↩︎

23. See childRmd/_11descriptives.Rmd file for other codes↩︎

24. See childRmd/_12crossTables.Rmd file for other codes↩︎

25. See childRmd/_18survival.Rmd file for other codes, and childRmd/_19shinySurvival.Rmd for shiny application↩︎

26. JAMA retraction after miscoding – new Finalfit function to check recoding↩︎

27. See childRmd/_23footer.Rmd file for other codes↩︎

28. Smith AM, Katz DS, Niemeyer KE, FORCE11 Software Citation Working Group. (2016) Software Citation Principles. PeerJ Computer Science 2:e86. DOI: 10.7717/peerj-cs.86 https://www.force11.org/software-citation-principles↩︎

 

---

A work by Serdar Balci

drserdarbalci@gmail.com